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FLOW OF MAGNETIZABLE PARTICLES IN TURBULENT AIR STREAMS 

by 

Kent Ritter Davey 

Sulanitted to the Department of Electrical Engineering on July 27, 1979 in partial 
fulfillment of the requirements for the Doctors of Science. 

ABSTRACT 

The requirement of particulate removal from turbulent flows arises in coal 
desulfurization, mineral beneficiation, water purification, particle research 
\diere particle loss is undesirable, and aerodynamic drag reduction, where the 
containment of particulate in a quasi-stationary manner within a turbulent 
boundary layer is desired followed by precipitation after a given length. 

Special consideration v/ill be given to particle precipitation. Light parti- 
cle (diameter < 2 Op) and heavy particle models (diameter > 2 Op) are developed. 

The first involves the numerical solution of a diffusion equation in which 
boundary conditions are imposed only ^ere particles enter the volume of interest. 
Inertial effects are unimportant. The second model involves a momentimi balance 
in Lagrangian coordinates augmented by a diffusion force. The diffusion term is 
a^ed so that the theory is consistent with inertia and diffusion dominated 
limits, and accounts for the effects of turbulent eddies in spreading particles 
in flight. Both models lump the effects of turbulent eddies into a measurable 
diffusivity constant. 

Precipitation experiments with light and heavy iron powders acted on by a 
stationary permanent magnet structure are correlated with the numerical model 
predictions. A useful degree of accuracy in predicting particle precipitation, 
as con 5 >ared to existing analytical models, is demonstrated for light and heavy 
particles . 

Positive correlations with data encouraged the use of the heavy particle 
analysis in examining particle flight in an aerodynamic boundary layer over a 
flat plate. The model predicted that with conventional permanent magnets, 90% 
of the injected particulate is collected in 5 meters and 4% would be lost. An 
alternative is described in which a travelling wave structure is used to contain 
and shuttle particles along in the boundary la;'er. In this mode, in vdiich 
particles continually interact with the wall, appears that particle loss can 
be greatly reduced. 
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I. Introduction 

The flow of magnetizable particles in a tuxbulent air stream in the 
presence of an in^osed magnetic field is of interest to NASA (National Aero- 
nautics and Space Administration) in utilizing the e:q)erimentally observed 
phenomenon of drag reduction produced by the introduction of particles in a 
turbulent boundary layer. The earliest observation of this effect were in 
the 1940 's [1, 2 , 3] and have recently been expounded on by Landahl [4, 5], 
Boothroyd, and Rosetti, and Pfeffer [6, 7]. According to Landahl, the particles 
dissipate energy in the small scale boundary layer eddies. Long, thin parti- 
cles are most effective in reducing drag. Other investigators believe that 
stabilization of small scale motion leads to a reduction in turbulent stresses 
near the wall, and to an associated thickening of the wall layer. In gas- 
solid suspensions, 10-60 micron particles yield significant drag reduction 
with a maximum reduction observed using 30 micron particles. Particles larger 
than 100 micron however, increase the drag [7]. NASA is interested in whether 
particles can be introduced into a tuibulent boundary layer, ducted along the 
skin of a fuselage, and then precipitated. Figure 1-1 (a) illustrates an 
hypothesized configuration over an airplane wing. 

The interest in particle convection migration is widely ranging. Soo 
[8, 9] has been working with multiphase flows for some time. The applications 
of and interest in such flows are quite varied. Some exan^les are as follows; 

(1) Determination of the proximity of the particle flow to stream 
motion in order to use solid particles as tracers in studying 
the flow of fluids [10]. 

(2) Determination of the diffusivity of the particles with respect 
to the continuous phase [11] . 
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(3) Determination of the diffusion of fuel and air as influenced by 
the relative motion of the fuel particles and air stream in order 
to optimize the combustion of solid field particles [12] . 

C4) Determination of the relation of the motion of solid particles and 
the field stream in sedimentation studies and pneumatic conveyance 
applications [13, 14]. 

More germaine to this investigation is work being done in the area of 
magnetic precipitation [15, 16, 17]. Liu and Lin give an excellent overview 
of work in this area. It is now clear that the use of magnetic fields in 
removing pyrite (desulfurization) and other inorganic contaminants from coal 
will be important in the next few decades. An IEEE Magnetics conference this 
past summer (1978) [18] , as well as in 1975 [19] , summarized the current work 
in this area. Currently, most magnetic precipitators use steel wool or a 
similar magnetizable mesh to enhance the magnetic gradient and precipitate 
particles (Fig. l-l(b). The first major contribution to this area where the 
effect of particle inertia was considered appeared in a thesis by Clarkson 
[ 22 ]. 

The central theme of this thesis is the flow of particles in turbulent 
air streams, and particularly in the boundary layer interaction. Towards this 
end, two practical avenues of research emerge. A detailed study of the precipi- 
tation of particles from an aerodynamic air stream is considered. Precipitation 
is required for removing particles from the boundary layer. It is also a 
measure of the turbulent diffusion magnetic migration processes at work in the 
flow. This thesis will focus on precipitation because of both its basic 
inplications and its practical application to drag reduction and particle 
management . 

(1) Small particle theory (< 10 micron) 
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Inertial effects are small but turbulent diffusion must be 
considered. The determination of particulate density is acconplished 
by solving a two-dimensional diffusion equation. 

(2) Large particle theory (> 40 micron) 

Turbulent diffusion is less dominant, and inertial forces are 
significant. Particulate distribution is analyzed through a 
numerical integration of the momentum equation, i.e. of determining 
the particle trajectories modified by the inclusion of a diffusion- 
type force. 

The principle contribution of this study is the incorporation of 
turbulent diffusion theory with an imposed magnetic migration process both 
with and without inertia effects. In this work, the imposed migration is 
magnetostatic; the nature of this magnetic force to particle interaction will 
be explored in depth. 

This thesis begins with a review of the pertinent background information 
on turbulent flows and the prediction of particulate diffusion. The nature 
of the particle magnetic force will be discussed and the inherent difference 
between electric and magnetic precipitation considered. Small and large 
particle concentration theories will be developed and followed by a presen- 
tation of the experimental apparatus, procedures, and theory correlation. 
Finally, the computer model will be used to simulate the flight of particulate 
in an air stream over a flat plate. The objective in this final study will 
be to understand the controllability of particle confinement in the boundary 
layer by altering the density profile. A steady density in the layer is 
desirable to obtain the drag reduction benefits. A brief consideration is 
given in the concluding chapter to using a traveling wave structure for 




Figure 1-2 Experimental Apparatus 
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walking particles along a wall in a boundaxy layer. Figure (1-2) shows sich 
a structure mounted on the t<^ of the duct. 

Before beginning turbulent diffusion theory review, a cursory examination 
of the basic experiment is in order. The apparatus used in studying both the 
precipitation and ducting of magnetizable particles is shown in Fig. 1-2. 
Particles are injected through a copper tube at A and blown down a six-foot 
duct. In the precipitation experiments, a permanent magnet structure serves 
as the source of a periodic static magnetic field source. The magnets along 
the bottom of the duct enhance the precipitation of magnetizable particulate 
on the lov;er plate. Particles not collected by this field leave the duct at 
D. A velocity profile grid at B is used to promote and control turbulent 
air flow. 
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II. MAGNETIC THEORY BACKGROUND 

A. Magnetic Force 

The first question that must be considered whether considering 
particle precipitation or ducting is "what Is the nature of the particle 
magnetic force?" A ferromagnetic particle of radius a, permeability w. 

In a field of Intensity If, experiences a total force 

T - 4 IT a^( — - 1 ) ^ V(lfr.ff) (2-1 ) 

2 + ^ ^0 

The reader is referred to Appendix A for a derivation of the Eqn. (2-1). 

The force expression (2-2) Is subject to the following two restrictions: 

1) The particle is much smaller than the characteristic length 
over which the field changes, i.e., the assumption of constant 
H" external over the particle's dimension is valid. 

2) Particle-Particle interactions are small. 

For the largest particles the author will be using (100 urn in a 
field structure wavelength of 5 cm, assumption (1) is quite valid. 
The second assumption is questionable when heavy precipitation 
occurs. The particles agglomerate into nair-like structures 
and enhancement of local field gradients undoubtedly occurs. This 
effect is only significant in the vicinity of precipitated ag- 
glomerates. (The author attempted to rest, let the amount of pre- 
cipitant to low levels.) 

The two field structures used in this work are shown in Fig. 2-1. The 
sinusoidal wave structure used in the ducting experiments excites a trav- 
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Figure 2-1 Field Structures Used 
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eling wave* the speed of which Is determined by the frequency and winding 
pole pitch. The field above the motor Is Laplaclan, decaying exponentially 
In the y direction with the same wave number k as the drive current T. 

The permanent magnetic square wave structure has an Infinite number 
of odd harmonics, again Laplaclan and falling off In magnitude In the y 
direction as the Inverse of the harmonic. Thus, a square wave field struc- 
ture with surface fle’d would yield In the upper half plane 


H 


Ho 2: 

m»l 

odd 


A. g-mky 
miT 


sin mkx 1. 


cos mkx 1. 


( 2 - 2 ) 


The field above such a structure is in reality not a square wave. Figure 
(2-2) shows the typical surface normal magnetic force density 1/4" above the 
permanent magnet structure. This field can be decomposed into Its Fourier 
components and represented as 


IT = I 

m 

m=l 


sin mkx 1, 


cos mkx i, 


(2-3) 


Using the first harmonic sinusoidal field of Fig. (2*l(a)) along with 
eqn. (2-1), one arrives at the particle force 


or 


F = -o' (2k) 


F = -o i. 


(2-4) 



Figure (2-Z) Normal Surface Force Density 1/4" Above Permanent Magnet Structure 
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where a - (4 it a^)( — - 1) Uq 

Note, the net force Is y directed only and Independent of x. 

For the more general Fourier expanded field of (2-3) F*F becomes 

00 00 

r Z cos(m-Z)kx (2-5) 

m®l t*l 

The force is obtained by taking the gradient of (2-5) and multiplying by 
the constant a' of (2-3) 


F = -a' 5T 
m=l 


00 

^ H H- 

m ic. 

Jl=l 


(m+t) k cos(m-5,) kx i 

(2-6) 

(m-)l) k sin(m-t) kx ? 

N t\ 


Equation (2-6) contains many components that will contribute little to a 
particle's motion. When m = Jl, all the x directed force components vanish. 
Furthermore, if m il, both the x and y force components have a sinusoidal 

X dependence which averages out for interaction lengths greater than one 
wavelength. This averaging is more effective when l and m are quite dif- 
ferent. Thus the significant force contributions occur when Jl = m and is 
y directed. 

i hI (2n,k)/2"«ky ? 

m ' 'e y 
m=l 


F = -o' 


(2-7) 
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The reader may wonder if there are other significant mechanisms for 
producing a magnetic force not accounted for by this model. Remembering that 
the force per unit volume is Uq?T*vJT, it is evident that the existence of 
a permanent moment m yields a force other than in the direction Non- 

col inear alignment of in and IT results in a torque on the particle and a con- 
sequent particle spin. Hysteresis would have the effect of giving rise to a 
non-col inear magnetization in a changing external field. Increasing tem- 
perature can cause permeability to decrease (true of most ferromagnetic 
materials) or to increase (e.g. magnetite below the Curie temperature). 

None of these effects will be considered in this thesis. The author will 
be using ferromagnetic powder (p - 7 x 10 kg/m ) with a permeability much 
greater than 

B. General Precipitation Remarks 

The author shall for the remainder of the chapter limit the discussion 
to particle flows where inertia is unimportant. The aim, specifically, in 
this section, is to gain an understanding of the important general proper- 
ties common to electric and magnetic precipitators. The added complication 
of inertia effects will be considered in Chapter 4. 

We shall consider the flow of particulate in a fluid of velocity u acted 
on by a migration force F. The particle flux is 

T = n(u + F/B) (2-8) 

where n = particle density 

3 = Stokes drag coefficient 6irnr 

and q¥ = electric precipitation 

F = 


a'7H«H = magnetic precipitation 
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Slnce 1^ + 7*T ■ 0, the Lagranglan density follows (assuming 7»u ■ 0) 
along (2-9) 


U + F/B 

When inertia is negligible, Eqn. (2-9) describes the state of the par- 
ticulate. The problem with (2-9) is that the velocity TJ in turbulent flows 
is a highly fluctuating quantity. A typical electric precipitation with 
negligible self field effect has a divergence free force field (V»qr = 0). 

In such an electric precipitator, the particulate density is constant along 
any trajectory flow line. One might wonder if the density would ever de- 
crease in a channel even when precipitation occurs. Equation (2-9) shows 
that whenever a packet of particulate enters the volume of interest, the 
density in the packet remains unchanged no matter how random or fluctuating 
its path. The particulate associated with a trajectory is removed from the 
flow when it meets a surface. Trajectories entering through solid surfaces 
(e.g. side walls in a channel) enter the volume of interest with zero density. 
It is the continued mixing of these zero density trajectories with initial 
particulate trajectories (which are themselves constantly being removed by 
precipitation) in turbulent flow that leads to a density decay down a channel 
with a constant cross-sectional value. This is in fact the Deutsch model we 
shall examine in the next section. The point is that all analyses can be 
said to be specializations of eqn. (2-9), the key issue being how to 
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handle the IT nf eqn. (2-9, b). Before applying this general theory to 
electric and magnetic precipitation in a channel, the author wishes to 
point out a general property of all magnetic systems. 

It follows from (2-9, (a)), that it is impossible to act on particulate 
with a force field whose divergence is positive definite, and have the dens- 
ity increase with time. This explains self charge spreading in an electric 
system where V*F is proportional to the charge squared. The author will 
now prove that the density always decays along trajectories for magnetic sys 
terns in curl -free regions ( in quasi statics this means no current density T) 

I 

The proof follows by demonstrating the divergence of the force field 
(**V(H»F)) is positive definite. A simplification of the force term can be 
made 


7(iT»H) = 2H*VH + 2li X (V xH) = 2H*VH 


(2-10) 


Using the Einstein summation convention, the term of interest can be written 


V*V(H.H) = 2 gj7 (H,. 




o 9 9 

9Xj 9x^ 


(H.Hj) 


( 2 - 11 ) 


The last step follows after applying V-H = 0. Applying the divergence free 
requirement to the last term in (2-44) again leads to 


9H„. 9H 

w * 


(2-12) 


9H. 9H. 

Now since V x H = 0, ^ and the proof of a decaying particulate 


density is accomplished. 
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C. Electric Versus Magnetic Precipitation 

Before examlnln In detail the full precipitation problem, it Is helpful 
to consider the nature of the.prec1p1tat1on process, and to Identify the In- 
herent difference between electric and magnetic precipitators. Electric pre- 
cipitators have been around for some time and analyzed extensively, but there 
lies 1n the physics of the magnetic precipitator a fundamental dissimilarity 
that warrants special attention. With this goal in mind, the author will 
compare the performance of a charged particle precipitator in the limit of 
complete mixing (the Deutsch model) and the limit of laminar flow with the 
magnetic precipitator in the same flow regimes. 

Electric Precipitation ; Figure (2-3) is the basis for modeling electric 
precipitation in a fully-mixed turbulent duct flow. Particles of charge q 
in a plug flow velocity U are acted on by a constant vertical field E^. It 
is assumed for the moment, that gravitational forces are negligible. The 
particles are charged before entering the precipitation region, perhaps by a 
corona source. The axial dependence of density is derived from mass conserva- 
tion arguments. In this mixing model, also referred to as the Deutsch model 
after its originator, the density is uniform over any cross-section because 
of the turbulent mixing. Balancing the particle flux entering with the flux 
precipitated and the flux leaving gives 

U(n(x + Ax) - n^x)) 2aw = — ^ w A x (2-13) 

or 

.i!a 

dn _ B 

3x " 2aU 


(2-14) 




The density then decays exponentially as 



exp - 


(X ■ Xq ) 


(2-15) 


where 



= 2a U 


The precipitation length p<j represents that length when roughly 2/3 of the 
Initial material has been precipitated. 

The equivalent precipitator in the laminar flow model is shown in fig. 
(2-4). An analysis based on mass conservation again follows but the trajec- 
tory flow lines are of importance now. The flux at any position in the duct 
is given by 




(2-16) 


Since the divergence of flux must equal the negative time rate of 
change of the density, it follows that 




- - 

’x- — 


(2-17) 


Note that the divergence of the T field should appear on the RHS of (2-17), 
ignoring this term is equivalent to assuming self charge effects are neg- 
ligible and that migration is dominated by the imposed field E. Any di- 
vergence-free force field will result in a constant density along the 
trajectory lines. 
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In the particle's reference frame, (2-17) can be written 


^ - 0 along 


dt 




(2-18) 


Equation (2-18) reveals that all particles are removed In the time It takes 
for particles to fall the height of the duct (t * The turbulent 

model never removes all the precipitant. The equivalent precipitation length 
In the laminar model Is 


2a(l-i) U 

S' qE„/6 


(2-19) 


Not only does the Deutsch model predict that the precipitant is never com- 
pletely removed, it predicts a longer precipitation length, their ratio being 



* 1 


i 

e 


( 2 - 20 ) 


The fact that these two lengths are comparable should be quite surpris- 
ing. In the Deutsch model, particles are supplied to the wall via the tur- 
bulent mixing and then precipitated because of the electric field. In the 
laminar model, migration to the wall is dominated by the electric field 
force. The precipitation lengths will provide a basis of comparison in 
examining the magnetic precipitation equivalents. 

Magnetic Precipitation ; The fully mixed magnetic precipitation model 
is shown in Fig. (2-5). Magnetizable particles of permeability y, radius a. 







A flux conservation 
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and density n enter the duct with plug velocity U. 
balance requires 

U 2a (n(x+Ax) - n(x))w ■ A x w (2-21) 

where F Is the magnetic force (Ignoring gravity). If the magnetic field 
structure Is a linear travelling wave, (2-5) gives the force as -ake"^*^^ 1y 
and (2-21) becomes) 


i dn . ake*^^^ 

n 3x ” “ UB(2a) 


( 2 - 22 ) 


Since the precipitation occurs at the wall where y ® 0, 


n . 


= exp - 


(x-Xq) 

P3 


where 


(2-23) 



2a U 
ak/B 


Again particles are supplied to the wall by the turbulent eddy mixing and 
precipitated via the magnetic force at the wall. Equation (3-23) reveals 
that the most effective precipitation occurs when k Is very large, i.e., 
when the traveling wave structure has a very small wavelength. The physics 
duscussed earlier would dictate that a small wavelength field would have a 
high gradient and thus a large particle force. 

In the laminar flow model (Fig. 2-6), particles of permeability y and 
initial density n^ enter the duct in plug flow with velocity IJ. A traveling 
wave source again exerts a vertical force ake’^*^^ 1^. The particle flux is 




Figure 2-6 Laminar Flow - Magnetic Precipitation Flow Trajectories 
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r 


U n V- 


0 



(2-24) 


Equating the flux divergence with the negative time rate of change of density 
gives 




^x- 


0 


Vn * - 


2ak^e*^*^^ 


(2-25) 


This Is a case where the Imposed force field has a limited divergence (see 
eqn. 2-9, (a)), unlike our electric precipitator. Thus eqn (2-22) can be 
written ] 

dn _ 2ak^ a’2ky 
dt ” “ S 


along 


dr ^ 
dt “ 



oy< -2ky j 

8 ® V 


(2-26) 


The decay is consistent with the spreading of the particle trajectory lines. 
By contrast in the electric case, the two-dimensional trajectories are paral- 
lel. Integration of (2-26(b)) gives the particle trajectory x and y depend- 
ence on time as 

X * Ut (2-27) 

y = ^ - M^) 


(2-28) 
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where the particle is assumed to begin its trajectory at x ■ 0 and y * y^. 
The density is found by substituting for y in equation (2-26(a)} and integ 
rating to give 




(2-29) 


The examination of these results as a function of starting position y^ 
is found by substituting (2-27) into (2-28) and setting y * 0 to give 


X 




(2-30) 


From this result or eqn (2-26, b) it is evident that there exists an op- 

dx 

timum wave number for minimizing precipitation length. Setting ^ = 0 in 
(2-30) gives 

2ky 

e ° (ky^ - 1) + 1 = 0 (2-31) 


The solution to this equation occurs when ky^ » 0.8. Because y^ can be no 
longer than 2a, optimum precipitation occurs for the wavelength 

A = ^ =1 m for a 12.7 cm duct. (2-32) 

One can now compare the precipitation lengths for the two limiting 
regimes using this optimum k as was done in the electric case study. The 
result is that the ratio of the laminar to Deutsch length is .55. This com- 
parison would give the laminar case an unfair advantage. The author wishes 
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to stress the difference of the above case study results from the electric 
analog. 

The electric case showed no difference (to an order of magnitude) of 
precipitation lengths for two entirely different flow regimes. The magnetic 
case gives entirely different results depending on the wavelength. The cor- 
rect precipitator design for a highly mixed turbulent flow where particulate 
Is supplied to the walls primarily by turbulent diffusion Is to Install a 
very short wavelength field structure at the precipitating surface. A good 
design In a laminar flow regime asks for a wavenumber roughly equal to the 
reciprocal half duct height. Too small a wavelength means little field pene- 
tration Into the duct volume, while too large a wavelength results In little 
field gradient and thus little force. This regime has no help from fluid 
motion to get particles to the wall. 

The foundation of two basic precipitation models is now laid. Further 
extensions to these models such as the Inclusion of additional field harmonics 
to more accurately represent the H- field could be made at this time, but the 
author wishes to Incorporate these Into the full precipitation models, A 
refinement of the magnetic precipitation models presented In this section is 
developed in Appendix B and the model predictions appear for comparison pur- 
poses with other results in Chapter 6. 
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III. TURBULENCE PHENOMENA 

The essential problem in effectively predicting the profiles of magnet- 
izable particles in turbulent air streams is to account for the effect of the 
turbulent diffusion. Three major camps have evolved in the study of particles 
in turbulent air stream: the flux conservation camp, the momentum or force 
balance analysis, and the stochastic or statistical attack. Much empirical 
work is intermixed with theory and as a result, many terms and seemingly, 
unrelated variables are used to represent effective turbulent diffusivities. 
The author wishes to introduce some congruity to these three avenues by an 
initial section defining the relevant parameters and discussing their inter- 
relation. 


A. Parameter Introduction and Interrelationships [1,2] 

In laminar flow, it is known that the following relations hold for 
viscous stress t and molecular diffusion D. 


T 3v 

P '' 3y 


(3-1) 


|^=07^n (3-2) 

where p = density of the fluid 

n = concentration of the particulate (1/m ) 
v^ * fluid velocity in the flow (horizontal) direction 
y = direction perpendicular to the duct walls 
V = kinematic viscosity 

By analogy, it is hypothesized that there exists a similar relation 


ship in turbulent flows, i.e.. 
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dV 

(v + Vf) -d^ “ V. 


3V. 


V ay t ay 


(3-3) 


II ■ (0 + 0^) V^n - D^V^n (3-4) 

Here v^, the turbulent or eddy viscosity (often given the symbol e) and 0^, 
the turbulent diffusivity (often given the symbol Sg) are much larger than 
the kinematic viscosity and molecular diffusivity respectively. 

The conventional starting point in most fluid mechanics derivations is 
to begin with the Navier-Stokes equations, split the velocity components into 
mean and fluctuating parts (e.g. v^ = 7 + v'), and then average these equa- 
tions in time [3]. The x component of the momentum equation becomes 


av 

X 

at 


_ av 


+ V 




av. 


X ax 


+ V 


y ay 


+ 1_ ( 7 ^) + i- {Tnrn = 1 + 

ax ^^x ' ay ''^x y-' p ax 


2 - 
vV U.. 


(3-5) 


where V»v is assumed zero, second order terms have been dropped, no z 

— 7 

dependence exists, and represents the average of the x directed perturba- 
tion velocity squared. In steady state this reduces to 

0 s - (v' V' ) + vV^v (3-6) 

ay X y' X ' ' 


The first term on the right hand side of (3-6) is hypothesized to be 
a 

ay ^^t 3y ' 


which completes the connection with equation (3-3). Since the ratio of 
is a constant again, by analogy, it is believed that D^/v^ would be a con- 
stant. Thus the ultimate goal of determining could be found if were 
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known. The major drive tjien is to determine v^. 

One of the earliest schemes for determining the so called Reynolds stress 
(v' V') was through the large scale eddy mixing length, given the symbol SL. 

A Jr 

The net flux of particles/area/time In the y direction assuming particles 
are transferred across a horizontal plane by the same perturbation velocity 
v; is 

Ty “ "(y + 4y) - n(y) • v^(n(y) + ^ |j-) - n(y)) 

or 

where the representative length In the y direction has been replaced by the 
Prandtl mixing length. The Prandtl mixing length Is effectively the mean 
free path of a pulse of liquid and Is thus a measure of the scale of the 
turbulent eddies. Comparing (3-7) with (3-3) and (3-4), one sees that 
y/yi = D|. if equals D^. Davies [2] points out that due to the elongation, 
both shear and angular, the processes of momentum and mass transfer will be 
effected. The mixing length i should be replaced by slightly larger length 
which would depend on the physical properties of the system and the intensity 
of the turbulence. A further extension of the mixing length concept is to 
relate it to the shear stress as follows: 

F = (1F>' 

giving 


(3-9) 
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Another parameter, the shear velocity U^, 



has been used to evaluate v^. 

(3-10) 


Here, Is the shear stress at the wall. The shear velocity can be cal- 
culated by determining the pressure drop between two points in a duct; the 
details are outlined In a thesis by Videla [3]. From relations (3-7), 
(3-8), and (3-10), It follows that 



wall 


wall 


Prandtl first suggested that mixing length was proportional to distance 
away from the wall through the Von Karman constant (k ~ 4). 

5, * ky (3-12) 

Michel et al . [5] and McDonald, et al . [6] have calculated relationships 
for mixing length as a function of distance from the wall, but their transi- 
tion across the boundary layer is questionable (see Fig. 3-1). Videla 
voices a caution about developing the mixing length approach. Apparently 
a more realistic transition can be framed by focusing on the character of 
the eddy viscosity and its associated shear stress rather than the mixing 
length concept. 

An approach the author feels may be more fruitful was first fostered 

by Taylor [7]. The argument runs as follows: the position y(t) of a given 

fluid particle is 
t 

y(t) = I v^(t')dt' 

0 


(3-13) 
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Viscous Sublayer - Laminar Flow Prevails 


Figure 3-1 Different Regions of Turbulent Duct Flow 
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where v^(t') is the Lagrangian perturbation particle velocity In the vertical 
direction. (We could just as well analyze the perturbation from the mean in 
the X direction.) Now 

t 

y(t) v^(t) * I" ^ y^(t) * j v^(t) dt' (3-14) 

0 

where the overbar signifies average (I.e., the Integrand Is first averaged over 
t). Finally, It follows that after a change of variables (T = t - t') and for 
long times 


1 3 

I at 


Tit) 



Vy(t) Vy(t-T7 




(3-15) 


2 

where v' is the turbulent intensity in the y direction obtained by squaring 
the perturbation velocity in that direction and then averaging on time. 

The symbol T is called the Lagrangian integral time scale and represents 
the time it would take for neighboring fluid eddies to become completely un- 
correlated. It is often tied in with the Lagrangian integral scale length 
and related to the integral time scale through the average turbulent velocity 
[ 8 ]. 


L = Tv^ 


(3-16) 


This length or Lagrangian eddy scale as it is also called, represents the 
length a fluid element would travel before originally neighboring fluid 
elements would become uncorrelated. It represents the size of the energy 
carrying eddies, those which get their energy from the main stream flow. 
The smaller eddies dissipate their energy through viscosity. Another form 
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of this length scale Is given as 
00 

L « I R(y) dy 
0 

where R(y) = velocity correlation = 


VyUj v;(y ) 



averaged on i 


The alternate correlation in space 


( 3 - 17 ) 






averaged on t' 


(3-18) 


Hinze [8] argues, must be equivalent to the space correlation when working 
with Eulerian velocities. 

Before drawing the connection between (3-15) and D^, it should be ob- 
served that in homogeneous turbulence (which is what the author will be 
assuming) the above analysis follows in all directions. Furthermore, the 
time and integral length scale will be the same. Finally, the above analysis 
was valid in a Lagrangian frame where tagged particles are followed. The 
measurements done in this paper will be Eulerian. Following a lead by 
Barfield, et al. [9], the author chooses to assume the Eulerian statistics 
are valid. The reader is referred to a thesis by Chadam for further dis- 
cussion of this topic [10]. 

The differential mass in the control volume of Fig. 3-2 is 


dm * p dx dy dz (3-19) 

Locally, the normal gradient of velocity is approximately a constant 
3v„/3y. The differential x-directed momentum is 
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Figure 3-2 Momentum Transfer In Turbulent Flow 
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3v^ 

dPx * p(“3]^) 


(3-20) 


Integrating over y gives the net x directed mor..dntum over this Infinitesimal 
rectangular area 
9v. 


'’x * 2 ^ ly" 


(3-21) 


Differentiating and averaging gives the net longitudinal stress as 
3v. 


_ 1 


'IP if 


xy 


(3-22) 


— 2 

Thus one observes from (3-15) that the turbulent Intensity v* and the La« 
granglan time scale are related to the turbulent kinematic viscosity. 



(3-23) 


One more point needs to be considered before the procedure Is consistent. 
Yermolll and Taggert [11] have shown after a lengthy study that the diffu- 
sivlty and eddy viscosity are not equal and indeed vary with particle size, 
density, and mixture concentration. A more accurate procedure would be to 
let 


= Sv^ (3-24) 

(see Videla [3]). However, their results showed 6 varied between .9 and 2 
for sundry sizes and mixtures. In addition, Chien and Einstein [12] 
showed that for various mixtures, 3 was almost nearly one. The author will 
in this work consider 3 = 1.0 henceforth. There is a fourth approach 
based on the time scale necessary to dissipate the energy of the eddies. 
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Jm-) but the Interested reader 1n referred to H1nze. 

V n 

To summarize, the different methods of approximating turbulent diffusion 
are as follows: 

(a) through mixing length 


2 ^''x 

Vi ■ ■ v' £ 

t 3y y 


(3-25) 


(b) through shear velocity U* 


Vi = U*£ 


(3-26) 


(c) through turbulent Intensity and the Lagranglan {- Eulerlan) time 
scale 


v^ * v^ T 


(3-27) 


All of the above assume that = $v^ = v^. It remains to demonstrate how 
these three attacks can be used to determine particle concentration. 

B. Flux Conservation 

The earliest work using this method for approaching turbulent diffusion 
was done by Schmidt (1925) and Prandtl (1926) [13]. A detailed analysis of 
the theory is outlined by Bauday [14]. We shall for the remainder of this 
chapter consider no other external forces aside from gravity. Beginning 
with the mass conservation equations, all quantities are broken up into 
steady state and fluctuating parts as was suggested earlier, yielding an 
equation of the form 

n' Vy + nWy = 0 

where y = 'vertical direction 
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v' ■ perturbation velocity In the y direction 
Wy ■ fall velocity in the y direction 

n » average concentration (3-28) 

n* » perturbation from average concentration 
Here Schmidt hypothesized the n* v' « ^ which yields the so-called 

Schmidt equation 

“t I7 * “y " ' ° 

The author wishes to outline a flux conservation argument because of 
its applicability to forthcoming work. If T represents the flux of particles, 
then the net particle flux T (part1cles/area/t1me) becomes 


f * convection + diffusion + migration (gravity) 


= vn + (D^ 


“molecular) 


(3-30,a) 


The divergence of flux must equal the negative time rate of change of concen- 
trati on. 


3t 


9(v„n) 3vln 3(v n) 

3x 3y 3z ' t 


+ D„) V^n 
m 


3{o)yn) 

1y 


(3-30, b) 


A simpler expression of (3-30) results from noting that the net y directed 
flux must be zero. Averaging (3-30(b)) and neglecting Dj^^olecular small, 
gives 


0t|+Wyn = 0 (3-31) 

The reader should note that three u'=sumptions are inherent in (3-31 )— steady 
state conditions exist, the particulate is dispersea evenly over the channel 
and i s the-zefore considered a continuum, and no net buildup of particulate 


occurs. 
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The work the author has examined 1s quite extensive and centers on mar l pu- 
latlons of equation (3-31) [15,16,17,18]. Graf and Raudkivi and Yalin al' 
give a broad scope of the general analysis, Csanady Is quite clear, and 
Hinze Is the most detailed and complete. Taggert and Yermolll as well as 
Videla are clear, supplying much of the background justification. 

The key to the problem Is In finding a diffusion coefficient to match 
a given turbulent flow. Assuming Is constant as Hurst (1929) [11] first 
did. 


n _ 


exp 


- [w (y-a)]/D^ 


(3-32) 


where n^ Is the concentration In the center of the duct at position y * a. 

A more realistic approximation to D^. may be gotten through (3-3). In 
a turbulent rectangular channel, the shear stress is [11] 


T = T (1 - Ji.) 
° ^0 


(3-33) 


where Is the stress at the wall and y^ is the distance to the center of 
the duct. Ippen [19] expressed the flow in a turbulent channel as 




(3-34) 
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Von Karman and Prandtl Independently determine the turbulent velocity rela- 
tion in 1934 [11] to be 
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(3-35) 
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using (3-34), (3-33) and (3-3), one can Integrate to get 


1 


w. 


JL . . yp-y ? kU, 

‘•yo-(Wk u*)-‘‘*y+v/k U* ■* 


(3-36) 


and if (3-35) is used. 


w 


n_ _ r’^o"'^ 3U*k 
‘•y«-b y* 


'b •'0 


(3-37) 


where n^^ is the concentration at some height b. 

Indeed the various attempts to obtain an accurate approximation for 
are quite involved. Videla chooses to evaluate different diffusion coeffi- 
cients in the different flow zones shown in Fig. 3-1. The notion that 
depends on distance from the wall seems reasonable because of the smaller 
scale eddies near the wall and the altered flow in the boundary layer. The 
summary of his work is as follows: 


Laminar sub-layer 
Buffer zone 
Wall zone 


Vt = V 


K u* y u* 

s ^{1 + JL - — tanh (a -rr-)} 


V 


,, = V y*- where i = non- 

t ^ dimensional 


vertical 
m * constant 


Defect zone 


Vt = U* y^ k y.( 






m' 


(3-38) 


Finally, White attempts to make correlations between and the Prandtl 
mixing length [20]. 

Several observations have sparked research along alternate lines [17]. 
For example, experiments indicate that the Von Karman constant (3-35) of 
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flow with suspensions of neutrally buoyant particles decreases with In- 
creasing concentration while the turbulent Intensity Increases. This runs 
contrary to the assumption that flows of neutrally buoyant particles should 
remain unaffected since no energy Is required to suspend them. 

C. Momentum Equation Attack 

The first major attempt to describe the motion of particles In a tur- 
bulent stream was by Tchen [21]. His starting point was the momentum equa- 
tion for the particle 


4 ^^3 ^ ^ 
T Pp 

1 


4 3 14 3 

6imr (vv ), + yirr + 2 T ’"'p (-Hr * It > 


_P' 




T 


.3 

P 


+ 6rl /irp^ 


t dv- dv„ 
i^t-t * 


dt' + Fe 


(3-39) 


"V- 
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The subscripts p and f refer to particle and fluid respectively. The terms 
are explained as follows: 

1) force required to accelerate particle 

2} Stokes viscous resistance force 

3) pressure gradient force in the fluid surrounding the particle 
caused by acceleration of the fluid 

4) force to accelerate the virtual added mass of the particle rela- 
tive to the ambient fluid 

5) force due to the fluctuating, non-steady state flow pattern 

6) external force field such as gravity or electromagnetic forces 
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The analysis Involves auto correlations of the fluid velocity, via a 
Fourier Integral switch Into the frequency domain. The result Is an expec- 
tation value of the square of the particle's position [22,8,17]. The theory 
Involves several assumptions, one of which requires that the neighboring 
fluid near the particle not change with time. This assumption reduces the 
Tchen analysis to a pedagogical exercise for turbulent flows. 

Khosla and Lederman [23] attempt to build on Tchen 's theory and alter 
the above limiting assumption. Instead of allowing time rates of change of 
the particle (gx-) to equal those of the fluid They hypothesized 

aip 


dv„ 

Ht^ ' 3t^ ^ Yv 


P 


(3-40) 


where y Is an empirically determined constant. The analysis again proceeds 
by going into the fourier frequency regime, and seeks to find the ratio of 
the diffusivities of particle and fluid. They conclude that the results are 
sensitive to temperature, density of fluid, and the presence of various par- 
ticles. These factors limit the method to low frequency and low speed tur- 
bulent flow. 

A momentum attack that has much credibility and links back to Taylor's 
early work [7] was first fostered by Chadam [10] and later elaborated by Bar- 
field, et al. [9]. With only the external force of gravity, the equation of 
motion of the particle is 

- 3 ^= a(t) - evp - (1 - ^) g (3-41) 

where 
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Pf 

^ s density fluid , 
density particle 

and a(t) ® the accelerative forces imparted to the particle by the fluid 
turbulence. If the particle was being dragged by the turbulent 
fluid, in a strokes type drag, one could represent this accelera- 
tion as 3v^(t), which is what the author will assume. 

Chadam integrates (3-41) directly but without including the effect of 
gravity. His analysis proceeds as follows: 


t 

Vp(t) = v^e"®^ + e‘^^ j e“®^ a(5)d? 

0 

“ T ■ i I 

0 

t 

+ jf a(5)ds 
0 

Squaring and averaging gives 

w 2 

(l-e-8‘) + i (-3 + 4e-S^ - 

r r It 

where A = 

t+At t+At 

I f a(?) a(e) d^ de 

t t 


and 


A = 


00 

» 

Vf(t) Vf(t-T) dT 

4 

0 


(3-42) 


(3-43) 


(3-44) 
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where the note after (3-41) has been used. 

With the analogy in (3-15) with Taylor's work, a connection to fol 
lows f 


1 d 
I3t 
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'large* 


y - 2?Tt 

X R- 

r 
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* /v^(t) v^(t-T)dx 

—R* 

= v'^ T 

(3-45) 
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Note, attention is fixed on the change in average squared position. The 
introduced in (3-42) is a turbulent perturbation velocity, the same as v' 
introduced earlier. The objective is to determine the fluctuation from the 
steady laminar type flow. 

The author would suggest another approach— to integrate (3-41). First, 
rewrite the basic equation as (g" * 0) 


dv' 

Tf " = 6v ■ 


(3-46) 


The homogeneous solution of (3-46) gives as the impulse re- 


sponse U.(t-T) of the system. With 6v' the drive, v' 

0 X|. 

t 

v’ (t) * f 8v' (t) dt 


becomes 


(3-47) 


Thus 


Xp=(t) 


I# 

g-3(t-x) gyi 


dt 


(3-48) 


0 -CO 
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interchanglng the order of Integration gives 
t 

Xp(t) - I 0 - v;; (t) dt 


(3-49) 


•00 


Barfield [9] (after erroneously writing (3-49)) claims that it is rea- 
sonable to assume that the turbulent diffusion process is Gaussian and has 


a variance equal to 
t 


I (1 - di 


and a probability distribution function 


(3-50) 


f(x„) » ( 


1 


) exp 


1 

2a2 


(3-51) 


P /(Atto^) 

Evaluation of (3-40) gives the same result as in (3-44) which agrees with 
Chadam. Barfield then generalizes this attack by including the gravitational 
field. Integrating (3-41) as before gives (from 3-49 but in the y direction) 


Vp(t) = 


f(l - [v'(t) - dT 


(3-52) 


or 


yp + ^(8t-l) = 


B 


V 

f [1 - V^(x) dt 


(3-53) 


where the lower bound on the particle's position is changed to zero for 
physical reasons. Barfield again argues that a normal Guassian distribution 
applies and the the probability distribution of the function 


G = y(t) + A (3t-l) is 
3 


7t2 


f(G) • (-7^) exp l- 




2o2 


(3-54) 


with a mean of 
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y---^(3t-l) (3-55) 

r 

and the same variance as In (3-51). 

The Gaussian assumption Is generally accepted especially where turbulent 
diffusion occurs In the atmosphere [18]. Indeed, fine particle plume con- 
centration analyses are developed from this view point. If one begins with 
the diffusion equation, the analysis Is as follows: 

1^ + ^ n s D^V^n (3-56) 

whose solution 1s 


n = 


n„ -xV4D*t 
e ^ 

2Ao^t 


(3-57) 


The spread variance is 


00 


- 1 
n. 


X n dx = 2 Da 


(3-58) 


This variance is Identical with Barfield and perhaps guided his reasoning 


along these lines. If one thinks of as xZ and remembering (in line with 

A P 

Id 2 

Taylor's work) that ~ plume spread from a point becomes 

X 

[8,18] (iin a uniform x directed flow) 


n. 


n = 




(x-v t)^ 2 2 

exp { ^ \ ^ 

H 2aJ Z,\ 


(3-59) 


where C. is the number of kilograms released per unit time at the origin. 

In line with (3-45) it follows that = v' t for short times (distances) 

A A 

2 

and is 2v' Tt for long times. The steady state source distribution is ob- 

X 


-63- 


9 “5 

telned by Integrating (3-59) from zero to Infinity with aj » 2v' Tt, 

A A 

O 9 

°y " 2Vy Tt* o^ ■ 2v^ Tt at large distances. Thus 
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The author wishes only to point out the connection to the motnentum 
Gaussian approach here. Chadam makes the connection between this probability 
distribution and the flux diffusivity model using statistical arguments (in 
the next section). However, one last point should be made concerning the 
connection between time and distance relationships of diffusivity as seen in 
sections B and C. Equation (3-44, b) points out that represents the cor- 
relation of the fluid velocity with itself, that is the time it takes a 
fluid element to forget its neighbor. The diffusivity is constant for long 
times when fluid velocities are uncorrelated. Davies shows that in the cen- 
ter of pipes diffusivity is constant and changes over the wall region, as 
shown by Fig. (3-3). Near the wall, turbulent eddies are more correlated 
than at the center. (This is not exactly the same for rectangular ducts.) 
This process is considered ergodic in the sense that time averages are re- 
lated to space averages. 'Herein lies the reliability of the formulas listed 
in section A showing space dependence. 

Lastly, a method for obtaining diffusivity in the main core of the 
channel based on the pipe friction factor has been fostered by Dhanak [27] 
and used in predicting the performance of electrostatic precipitation [28]. 
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Dhanak begins by expressing the eddy viscosity In terms of the Fanning 
friction factor f of the channel, where 



This factor Is determine by either pressure differences down the pipe or the 
Blaslus's equation. Although the formula derivation Is uncertain, Dhanak 
concludes that 

D. * .0708 R. v^T (3-62) 


He compares the above prediction with 


V V 

D, . (3-63) 

3y 


(Which agrees with (3-6)) and finds good agreement. 

D. Statistical Approach 

Chadam [10] closes the loop between the Gaussian normal approximation 
outline in section C and the diffusion equation. He defines a function 
W(vi t,v^,t^) to be the probability of a fluid element with velocity at 
time tg to have velocity v' at time t. A similar function W(x,t,XQ,tjj) is 
defined for transition of particles from position x^ to x. Csanady states 
that most turbulent flows behave according to a Markov process, one in 
which the velocity autocorrelation 


vi(T)vj^(T-t) 




R(t) = 


(3-64) 
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where the numerator overbar Indicates an average over t. 

Chadam shows that In such a process » which can be assumed Gaussian, 


W((v',t,v^,tQ) * ' Y 

2tt I (1 - 

where A Is defined In (3-44) 
Furthermore, 


a (1 . e-28t) 


(3-65) 


W(x,t,Xg,tQ) = 





(3-66) 


which would agree with Barfield et al.'s equation (3-52). His next step 
however Is to find the probability distribution P(x,t) of finding a particle 
at position, and time t using 


P(x,t) = 


P(x‘ ,t' ) W(x,t,x' ,t‘ ) dx' 
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It is shown that the function 


t+At 


B(At) = 


a(C) dC 


has a Gaussian expectation value 


W(B(At)) = — — yT2 
(47rAAT)‘'^ 


iimx 

4AAt 


(3-67) 


(3-68) 


(3-69) 


The following mathematical Taylor expansion 



P(x'.t') + (x,t‘) ■ P(x,f ) - /Ax(P'W + PW) d(Ax) 

+ j ^ (PW + 2W'P' + P"W) dAx 
+ 0(Ax)^ (3-70) 


(where the primes in the integrals indicate space derivatives) 
reduces with (3-67) to 



I h 

8 3X 


(3-71) 


Thus, again a turbulent diffusivity results equal to v' T. 

The author wishes for completeness to examine two more approaches. 
Batchelor [24] defines a characteristic functi on (f)(^) in a normally dis- 
tributed Gaussian field 


<}>(?) = exp j 4 

(adopting Einstein summation convention) 

Conditional probabilistic arguments are then used to link <j)(c) to 
P(r,t). This involves transforming between ^ space to x space. The details 
are given in Barfield [4]. The result after much calculation is 



and 
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R^j(i) + Rj^(T) dt 


(3-74) 


Vi(t)v.(t+T) 

where ^ (averaged on t) 

vpr 

Now for large times the Integral In (3-74) becomes the Lagrangian time 
scale T and (3-71) becomes 


= vj vjT V^P(r»t) + ^ • V P(r,t) 


(3-75) 


Finally, Soo [25] begins by expanding the fluid velocity in an infinite 
spectrum of harmonics (isotropic turbulence) 


v' * v' 

X > 
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(3-76) 


He then couples this into the force equation 


dv 


~St ’ ■ ''p' * 9 


(3-77) 


Next, (3-76) is squared and averaged to give 
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After expressing the correlation function R(t) In terms of these Infinite 

X 

components and the scale of turbulence / R(x) dx. the Einstein equation for 

0 

diffusivlty [26] Is used to relate the diffusivltles of one phase of the 
flow to the other phase. The analysis is Intended for Intimately mixed or 

multiphase systems. His conclusion Is that the correlation, scale, and In- 

\ 

tensity of one phase of the system can be calculated from statistical re- 
lations on the other phase. 

The author will henceforth adopt the position that the turbulence ef- 
fects can be Incorporated Into a diffusivlty model where ^oiffuslon “ 
and Is represented by the turbulent Intensity v' times the Lagrangian 

time scale T = / | V' , t , t-T) ,, dt ^ physically related to 

0 v'^ 

the approximate time for turbulent eddies to decay. It Is related to the 
integral scale L (size of largest eddies) through the mean flow 

T = 4: . (3-79) 

''x 

has been found to behave roughly as b + d In channels. Its measure- 
ment will be described in Chapter 5. 
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IV. TURBULENT PRECIPITATION THEORY 
A. Introduction 

TVro precipitation models have been pursued in Chapter 2 : 

(1) a fully mixed tuxbulent model (Deutsch) where the particulate is 
sillied to the wall via eddy mixing and the cross sectional density 
is uniform (representing an effective for transverse diffusion 
equal to infinity) 

(2) a trajectory model where the flow lines are determined from particle 
conservation and in this method of characteristics analysis the 
density is found to decrease along trajectory lines. 

This chapter outlines a method of predicting magnetizable particulate 
precipitation from a turbulent airstream. Two analyses emerge based on the 
above models. The first, appropriate to li^t particles, applies where 
particle inertia is negligible. The theory is analogous to the Deutsch model 
so far as the physics of the problem (both accent the inportance of turbulent 
diffusion), and the analysis parallels that in the trajectory model. The net 
particle flux has four conq)onents : 

gravitational magnetic 

r = convection + migration + migration + diffusion 
= ^ n + a' n + DVn 

(3 P 

where D - Molecular ^turbulent ' ^turbulent 

The diffusion equation results after setting the divergence of (4-1) 
equal to zero. 

The second heavy particle analysis considers particles in which inertial 
effects are inportant. In -the consideration of very heavy particles, the theory 
to the laminar flow study of chapter 2. The particles with large 
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inertia ride throu^ turbulent eddies, and diffusion can be neglected. This 
analytical theory is expressed by the particle momentum equation vihich pro- 
duces infoxmation of particle trajectories: 


md7 


dTTTir 
P P 


6irnr + mg + a'viK’H) 
P ^ 


(4-2) 


The iji^ortant point of this investigation is in connecting these two 
analyses to produce an explanation of particle precipitation. Summarized in 
Table (4-1) is the interconnection of particle size, fluid flow, and the two 
analyses. The heavy particle turbulent section is split to indicate the model 
dependency on the degree of turbulence (lower graph) . In reality there exists 
a gray area in which both diffusion and inertia are important ; for this , a 
hybrid model will be developed in Section D- 2 . 

A synopsis of this chapter is shown in Table 4-2. The decision of vMch 
size-dependent analysis is valid must proceed from consideration of the rele- 
vant system characteristic times. Three models evolve in the light particle 
diffusion analysis, llie three differ by diffusivity representation, boundary 
conditions, and ccmputer simulation used. The heavy particle analysis is 
siibdivided according to whether diffusive effects are considered. 

B. Characteristic Times 

The major question to be answered is, "For what size particles are the 
effects of diffusion and inertia inq)ortant?" This is best answered by 
examining the characteristic times of the system. In balancing the inertial 
and viscous teims in equation (4-2), one arrives at the inertial -viscous 
time 

\.659 X 10'^ sec iron particle 

= —HL = i = (4-3) 

^ in- vis " 6Trnr ” ^ 


.3137 sec llOy iron particle 
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Degree of Turbulence 
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Llght Particle Analysis 
Diffusion Effects 


Characteristic Times 



Heavy Particle Analysis 
Inertial Effects 



Inertia Model 

Momentum Equation Attack - particles 
ride through turbulent eddies-no 
diffusion (like laminar flow theory) 


Hybrid Model 

Hybrid Momentum Attack - diffusion 
term added to momentum equation 
by superposition 


TABLE 4-2 


Chapter 4 Synopsis 
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For processes involving times shorter than experiments show that 

inertial forces must be considered^ a system vdiere the characteristic times 
were longer tlian can be shown experimentally to be accurately 

represented by the viscous dominated flux expression (4-1) . 

Ihe choice of proper characteristic system time deserves careful con- 
sideration. Because the particle's residence time is defined to be how long 
it takes a particle to traverse half the height of the duct, special attention 
to force terms of the RHS of equation (4-2) is in order. The basis of ccmi- 
parison is the time it takes a particle to fall the duct half-height in a 
viscous dominated environment. The gravitational- viscous time is 


grav-vis 


a 

mg/3 


3.905 sec 8u particle 


.02066 sec llOp particle 


(4-4) 


The equivalent magnetic-viscous time calculation requires more information 
concerning the field strength and v/avelength. For a single harmonic sinu- 
soidal field, the magneto-viscous time becomes 


T 


P 


0 





1 ) 


(4-5) 


and with a 
gives 



•)(}| — 1)vi^Hq , a 1000 Gauss, 8 cm wavelength field 


T 


mag -vis 


^102.5 sec 8u 

■ 

.542 sec IlOy 


(4-6) 
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Finally, the analogous diffusion time is 



a a 

^ r 

n 


(4-7) 


The measurement of the turbulent coefficient (-D) discussed in Chapter 2 

is outlined in Chapter 5. A representative diffusivity for sudi a system is 
2 

.002 Thus, the diffusion time becomes 


Tjj » 2.016 sec, independent of particle size (4-8) 

Examination of the times in equations (4-5) to (4-8) allow a more in- 
telligent handling of equation (4-2). For the 8u particle, the inertial’ 
viscous time is much shorter than the characteristic migration or diffusion 
time, indicating that the viscous -dominant flux analysis posed in (4-1) is 
valid. In the heavy particle analysis, however, the characteristic inertial- 
viscous time is not short. Therefore, inertia has a significant influence. 
The fact that the diffusion time is much longer than either of the migration 
terms, however, illustrates the reduced effect turbulent diffusion has with 
larger particles. The momentum equation (4-2) without the diffusion term 
appears to be the correct approach here. Thus there exist two well defined 
avenues of analysis depending on the size of particle being used and the 
resulting relative time constants. 

The larger particles used in the e:q)eriments were 55-65 microns in 
diameter in which the gravitational viscous time and magneto-viscous time 
are .082 seconds and 2.08 seconds respectively. Clearly these particles fit 
into an intermediate domain where diffusion is not negligible , especially 
over distances smaller than the duct half iieight (e.g. the injection tube 
diameter) . This intermediate region will be considered in section D of this 
chapter. 
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C. Diffusion Model - 4y Particles 

The above analysis revealed that 4y particles are viscous dominated, 
inertial effects are negligible, and thus a flux, particle conservation approach 
using equation (4-1) is in order. Assuming steady state opeiation, the 
divergence of T can be set equal to zero to give 

-V* (D^) + n - ^ ° 

Using a single haimonic field with tlie magnetic force equal to ake (4-9) 
becomes 


n mnl_^ . ( koie mg ^ 3*1 . .. 3n . 2o(ke _ n 

-v.(Dvn) + ( g f-l * u g - 0 


2-2ky 


(4-10) 


The next step in the analysis is to check whether any siuplifications can be 
made. In this study and indeed for most boundary layer analyses, the con- 
vective horizontal flux dominates the corresponding diffusive flux. A 
characteristic density analysis elucidates this point. 

Given a unifom injection of particulate during a time sufficient for 
10 grams to be spread over 2 meters of duct, the average concentration is 


(.127)^2 

and assuming a linear axial 
by 


^ (4-11) 

concentration profile, the (iensity n is typified 


n - -0.6 (x-2) ^ (4-12) 

^ 2 

In a 4.3 ^ flow with a diffusivity of .004 ^ , the convective and 
diffusive fluxes are respectively 


r 


convective 


1.3-t^ 
(m -sec.) 


(4-13) 
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^diffusicai * “ *0024 — 




(4-14) 


(m -sec) 

Thus the axial diffusion tem is justifiably ignored. As shown in the last 
section » vertical diffusion must be considered to be at least con^arable to 
gravitational and magnetic migration. 

The equation to be solved then reduces to 

_ |D 3n _ m£ ^ _ oke’^^ 

7 av av ft av "(iT 

ay 

. 2ak^e'^. 


- a£2“ . S&jn . * U (a-ly-a| an 

^ sy ay 6 ay 6 ay maxV a / ax 
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(4-15) 


where D is a function of y and x. Including the field harmonics t (4-15) 
becomes 

1 

-d 2^. ”£|il- V fy) faiix^y |a . v' (y)n - 0 (4-16) 

^ ® W n»sDc ^ a / ax ^ 


where 


Vjj^(y) » vertical magnetic migration velocity 
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aVm _ 4 __3 i »u 1 N y 3 - y 2k 
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00 
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The iii^osition of the proper boundary conditions on equation (4-15) is 
necessary for numerical solution. The fact that (4-16) is first order in 
X and second order in y suggests two boundary conditions on y (one each at 
the top and bottom of the duct) and one on x (at the entrance to the duct) . 
The boundary condition on x is simply the specification of the incoming 
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particle density (n^) . The principle of conservation of particle flux must 
be enployed to arrive at the vertical boimdary conditions. 

Further analysis is best pursued in stages. This step progression is 
needed to highlight the internal issues involved in a self-consistent solution. 
First, a non-causal theory is developed assuming a uniform diffusivity across 
the duct vdth a step change to zero a distance A from the walls. The 
equations are then represented in difference form. Axial differentiation 
will be represented by forward one-sided derivatives yielding a very straight- 
forward solution. This is then followed by two more exact causal theories 
in vMch representations of diffusivity are used. Boundary conditions and 
representation of axial and vertical derivatives are examined carefully to 
be consistent with causality. 

1. Non-Causal Diffusion Model 

The measurements of turbulent diffusivity (Chapter 5) along with the 
observations of Davies discussed in chapter 3 indicate the diff’isivity is 
approximately constant over the duct cross-section and decreases to zero 
near the walls. If turbulent contributions to the particle flux at the 
walls are negligible, a resonable model is one with the problem divided into 

(1) a boundary layer region close to the wall where diffusion is 
insignificant. 

(2) a region over the duct interior where turbulent diffusivity is 
constant with y, but does has an axial dependence. 

Thus, in both regions, the ^ term in (4-16) drops out. 

The above assunptions (summarized in Fig. (4-1)) along with flux 
conservation suffice to give the vertical conditions for the inner diffusion 
region. At. the lower boundary (y = A) , flux continuity demands 
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3n 
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- mgn + 


a'7(H*H)n] 

y ■ 


[-mgn a’V IT*H3 

y - A. 


or ^ . 0 at y - A^ 

At the t^per boundary (y - 2a - A) , continuity again requires 


(4-17) 

(4-18) 


[-mgn + a'V IT-H) - [-D^ ^ - mgn + a’V IT*H) (4-19) 

y - (2a- A) ^ “y y - (2a- A). 

Equation (4-19) iiqplies 

a . 0 at y • C2a-i). (4-20) 

In the layer regions it is appropriate to apply flux conservation, but 
with » 0. This is identical with the particle trajectory theory discussed 
in Chapter 2. 

Each of the above model divisions will now be examined in more detail. 
The greatest difficulty is with the lower boundary. There are two physical 
conditions that must be realized in the lower region- -no particle reintrain- 
ment occurs, and diffusivity for iron particles in the laminar sublayer must 
be zero (i.e. molecular diffusivity and brownian motion do not enter for the 
size scale particles we are concemed with) . 

An understanding of approaches enployed by hydraulic and chemical 
engineers to the precipitation problem is helpful. A method used often in 
hydraulic's literature (Barfield, et. al. - see Bibliography, Ch. 3) for 
modeling an absorbing wall is to force the density to zero at the wall. The 
pertinent diffusivity is considered to remain high near the wall in sedi- 
mentation research, for instance. This zero density requirement inplies two 
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conditi(sis \diich give credence to this assigitrent. 

(1) zero density at the wall insures a positive gradient away from the 
wall and thus only downward diffusive flux. 

(2) large migration (gravitational) flux to the wall with abnoimal 
depletion of the core region is prohibited by the assignment n ■ 0 

at the wall. The density required to be small in the neighborhood 

of the wall, since any large gradient would constitiite a correspondingly 
large diffusive flux (which could not exist in steady state) . 

A typical engineering problem is the dissolving of a salt at a wall into 
a flowing liquid. The wall molecular diffusivity is very inportant again. 

The wall for such a problem is often siq)plying flux into the duct or keeping 
the particle density in equilibrium, in which case the diffusive and migration 
fluxes balance (iii?)lying a negative density gradient at the wall) . At the 
laminar sublayer where the diffusivity increases tremendously, the normal 
density must be nearly zero. 

Appendix D shows the analytical solution of a three region diffusion 
problem (Fig. 4-1) with only gravitational migration, and no x dependence. 

The boundary conditions are n * n^ at y ■ 2a and n = 0 at y = 0. The 
diffusivities are allowed to be nonzero in the laminar sublayers. The 
result is that the normal density gradient at the boundary layer equals the 
ratio of the diffusivity in the sublayer to the diffusivity in the core. 

Also plotted are density profiles for various ratios of core to layer 
diffusivity as well as variation effects of diffusive to gravitational 
migration times. The density gradient in the sublayer is determined from 
trajectory theory and is fixed by the gravitational and magnetic migration 
terras (Appendix B) . Negligible layer diffusivity therefore in^lies zero 
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cors gradient. 

The condition can be examined less rigidly through the exponential 
magnetic force decay in the duct. Changes occur less rapidly axially than 
transversely. Snail axial flux changes imply a somewhat constant transverse 
flue. CHiis hypothesis will be used to explain density profiles in diapter 
6 ) . In the core , magnetic migration is small . One mechanism by vhidi flux 
is conserved from the core to the layer region is by diffusive flux di- 
minishing at the lower layer where magnetic migration increases. The 
diminishing flux appears as a decreasing gradient near the lower layer. 

Figure (4-2) shows characteristic profiles of these two cases as well 
as a calculated density profile for magnetizable particulate deposition vdiere 
the wall is absorbing (i.e., no reintrainment occurs) and molecular diffusivity 
in the laminar sublayer is zero. Zero sublayer diffusivity is a reasonable 
assun 5 )tion for 2-8 micron iron particles. 

From this there appears to be a sound basis for approaching the problem 

as a multi- region one, the core being a constant diffusivity diffusion 

problem wi.th boundary conditions 5 — « 0 at y « A and y = 2a - A. Differences 

<jy 

with other engineering disciplines have been considered. 

There now remains the task of matching the core region solution to the 
laminar sublayer regions. Because D = 0 in the sublayers, the method of 
characteristics approach discussed in chapter 2 can be applied. With a single 
harmonic field the rertinent equation is 

dn _ 2 ak^e’^^« dr _ rr ^ mg ake”^^^^ r/. 

^ ^ n along 3 ^ - U + §- ^ C4-21) 

and with the full harmonic field (4-21) becomes 



(4-22 (a)) 



Figure (4-2) Density Profiles 
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dr 

3t 


U 


max 





(4- 22(b)) 


Here the e:q>licit foim of IT has been inserted. Integration of (4- 21, a) shows 
that along the particle trajectory lines (4-22 ,b), the density decays as 


V (v ) 
n , mg ♦ m^^o^ 

n mg Vjjjcyj 


(4-23) 


where n^^ and y^ are the initial density and position of the particulate in 
the layer (i.e. at y * A). The particulate distribution on the bottom of the 
duct is given by (4-23) with three specifications -- n^^ is set equal to the 
density of the diffusion analysis ((4-16) at y = A) , y is set equal to zero, 
and the appropriate axial position for each density is found by numerically 
integrating (4-22 ,b) across the boundary layer. The splicing of the different 
region solutions is thus accon^lished. 

The i 5 >per layer must be investigated from a different perspective, 
because all vertical migration is downward. Since the diffusivity goes to 
zero in the upper layer, no mechanism exists for particles to enter this 
region. In other words, the particle trajectory lines emanating from the 
ipper boundary must have zero as their initial density. The upper layer 
correct boundary conditions are : - n = 0 and ^ = 0 . Although numerical 
integration shows that the two conditions give nearly identical precipitation 
results, the densities do differ throughout the volume. The zero density 
condition is only a result of no upward migration forces . The zero gradient 
condition exists independent of migration force direction- -this latter 
condition therefore seems to be more fundamental. Figure (4-3) sxmniiarizes 
the system representation and boundary conditions used. 
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Before examinixig the numerical method of solution, it is helpful to 
normalize the defining equati(xx (4-16) . Variable normalizations and 
relevant times are as follows: 
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i' T 


Diff 




■^Diff. 
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The diffusion equation is now written 


3^ . '^Diff Kn) . '^Diff aCn) ^ ^Diff 3 (n) 

' ZJ ~ T 4y T .. 3y ax 


grav-vis 


mag-vis 


•^Diff 

^mag-k 


^Res 


n = 0 


C4 


The boundaiy conditions are that (i^ = 1 at x = 0; ® Z. “ ^ » 


and ( 


f MbI = 0 

9Z 


atZ=l 


or n = 0 


-24,a) 


-24 ,b) 



-89- 

Method of Solution 

The proper method of solution in the duct diffusion region with an 
initial boundary condition on x is to sequentially step in x space. A finite 
difference system was used to numerically solve equation (4-15). 

Figure (4-4) shows the partitioning of the diffusion region into a number 
of points separated vertically by Ay and horizontally by Ax. In finite 
difference fom» equation (4-24) at point ij becomes » 


piin,j * a.i-i,j 


T, f 1 

' A 2 

J 

Diff \ 'Tgyav-vis 




V ^ . 

' '^Rss 

i dx 1 


T . i 

mag- vis ' 

“^Diff „ 

T — i 1 

^mag-k 


0 (4-25) 


Because all the first column elements are known (h , = 1), expressing 
(4-25) for the (j-2) elements yields a solution of the (j-'O elements. The 
net system solution is obtained by stepping column-wise from left to right. 
At the lower boundary if one were to specify that the next row of elements 
(i+3) have the same value as the (i+1) row elements of the same column, then 
the condition ^ ’ 0 at ^ would be insured. Matching the (i-2) row 
elements \d.th the (i) row elements at the upper boundary insures the same 
condition on the normal density gradient at the top. If the alternate 
condition (n*0) is desired at y = 1 - th«i the last set of difference 
equations must stop within the diffusion region, i.e., at row (i) . Appendix 
E lists the programs KPDI4 used with n=0 on the ipper layer, and KDDI9 with 
^ = 0 at the ipper layer. The correlation of the experimental result with 
the above predictions is discussed in chapter 6. 

A point about stability is necessary before concluding this approacli. 
The worst possible solution perturbation is oscillatory in y having the form 
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Figure 4-4 Numerical Solution of Diffusion Region 
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at a given grid point (see Fig. (4-5)) 


The criteria for stability is that such a perturbation not grow in x. With 
such an oscillatory pattern, equation (4-25) including diffusion and con- 
vection tenns only, demands that 




+ f . 


b.t 



^Diff 

^Res 


Ax 


(4-27) 


or vdth equation (4-26) 






^Res 
■^Diff ^ 


(4-28) 


If the pertuibation f is not to grow the second tern in brackets must not 
exceed 2.0 in magnitude. This condition translates to a restriction on 


i.e. 


Ax < 


^ ^Diff 
(2)*'^Res 


(4- 29) 


For the present analysis with a typical residence time of 0.23 secs, a diffusion 
time of 4.03 secs (D^ = .004), and Ay » .1, Ax must be less than .087. In 
this analysis. Ax was set to .0125. 

It is interesting that the first solution of (4-25) involved finite 
difference equations as above, but with two boundary conditions on x. A 
boundary condition of n = 0 was inposed approximately 24 duct widths down- 
stream. The finite difference system was in^osed as above , but the network 
of equations was solved simultaneously using a Gauss -Jordan elimination 
technique. Except in the neighborhood of the downstream point, the solution 
was of the same order as the above stepping procedure. The stepping procedure 
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eliminated the need to solve simultaneous equations but the question of 
causality when using the one-sided foxward derivative in x must be considered. 
The results of the non- causal deposition prediction are compared in Chapter 
2-A. 

2. Causal Perfunctory Model 

Three modifications of the non-causal solution are in order. The 
question of causality points to the first refinement, calling for a more 
careful examination of the boundary condition imposition as well as repre- 
sentation of first order partial derivatives. With no diffusion, equation 
(4-24,b) becomes a first order partial differential equation which can be 
solved using diaracteristic trajectory theory. The correct analysis in this 
limit would impose one boundary condition wherever particle trajectories 
entered the region of interest. Spatial derivatives in x and y would be 
evaluated in a one-sided difference representation always using particle 
location information of earlier time. The exact solution should be con- 
sistent with the trajectory analysis and in fact, degenerate to it as 
diffusivity decreases . One approach to the problem then would be to extend 
a finite difference grid through the i;?3per and lower laminar sub -layers, 
iiiQsosing the boundary condition n=0 at the iq)per surface. It would also be 
in order to use one-sided partial derivatives consistent with causality for 
non- diffusion terms. This solution would surely be more fundamental. 
IMfortunately, the refinement of the grid necessary to develop infoimation 
throu^ the boundary layer presented numerical difficulty to the small PDP 
11-03 coii5>uter available in this study. It was decided to continue to split 
the problem into regions. 

The nature and position of the split leads to the second solution refine - 
ment--the transverse dependence of the turbulent diffusivity. Equation (3-38) 
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reveals one representation of the diffusivity vertical dependence. The 


point was made in chapter 3 that according to Davies, the diffusivity is well 
represented by Figure (4-6) . The turbulent diffusivity is linear over 1/4 
of the duct half - hei^t dropping very quickly to zero in the nei^borhood of 
the laminar sublayer A. The thidkness of the layer is confuted from Sdilichting 
(chapter 5-[l]) to be 

A • .031 cm (4-30) 

y I 

V. 03325 U ^v^a ^ 

>1 avg 

Here U^yg» the average duct velocity, was assumed to be 4.57 ra/sec. Sdilichting 

also shews that the axial flof in the sib layer can be represented as 

7 1-1 

(.03325 a 

U- 22S (4-31) 

V 

The equation to b e solved in the core is identical to (4- 24) , b ut with 

the term of (4- 16) 

2 

^ 2L ‘^Diff ^ ^Diff M 1 + “^Diff an 

■ ^ ' Wvis Wvis FoiSJL 9=^ 

- l2iS_n • 0 (4-32) 

^mag-k 

The normalization of (4- 24 b) ■^as used in (4-32). 

The model alterations incorporating these refinements appear in Figure 
(4- 7) . Boundary ccxiditions are inposed on the top and left side of the duct 
where trajectories would enter. The grid is extended into the upper layer 
where the condition n=0 is inposed. 

The lor er boundary condition at y=A is necessary if is non zero at 
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that point (2nd order in y). As figure (4-6) points out, a discontinuity in 
exists at y*A. The imposition of D^-0 at y»A ijiplies significant alter- 
ations that will be discussed in the following section. The discussion in the 
causal model section of the boundary condition * 0 at y«A applies in this 
model as long as is non zero throughout core. 

Hie trajectory model was then meshed with the diffusion solution at y=A 
as explained in the approximate model solution. 

The finite difference representation of (4-32) consistent with causality 
and degenerating to the trajectory analysis as goes to zero becomes 



^Diff 


n. . = 0 


j-k 


(4- 33) 


Note the ~ term uses a two-sided derivative of density, while the migration 
and convection terms use a one-sided derivative. The solution is quite 


general in that it is applicable to any type of migration. If migration were 
up in some region and down in others (e.g. due to positioning magnets at 


y=2a as well as y=0) a point by point check of the direction of local migration 
would be necessary before one would know whether to choose n. , . or n. . . 
for evaluating 

9X; 

Causal representation of the x derivative convective teim necessitates 
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tlie sunultaneous solution of all densities in a single column before proceeding 
to the next column. The point to remember is that the axis in which dif- 
fusion is Important will require simultaneous solution of all points along 
that axis. The other axes give the directions of sequential numerical 
stepping. 

Appendix F lists the program KDPFZ used to predict the causal perfunctory 
li^t particle diffusion theory solution. The above discussion illucidatos 
the requirements causality in^osed on boundary conditions and numerical 
differentiation. A conparison of this models ' precipitation prediction with 
experimental data is in chapter 6-2. 

3. Causal Fundamental Diffusion Model 

As shown in Fig. (4-6), there exists a discontinuity in diffusivity at 
y«A. If is set to zero at the lower layer, the core analysis requires 
no lower boundary condition. The diffusion equation representation for the 
final row of elements in Fig. (4-7) involves only elements above and to the 
left of their location (i.e., backward in time). 

Except for the change in and this boundary condition, the analysis 
is identical to the previous section. The program listing appears in 
Appendix G. The theory precipitation correlations with data follow in 
Chapter 6A. A conparison of these two causal solutions differing by the 
lower boundary condition is also discussed in Chapter 6-A. The zero boundary 
condition theory predicts a 5-10% lower particulate deposition than the zero 
gradient causal theory. The most surpilsing result was that the no lower 
boundary condition causal solution preda ted a density profile identical in 
shape to the zero gradient theory, including a zero density gradient at y=Al 
The conclusion is that this third no lower boundary condition causal solution 
is more fundamental. 
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D, Heavy Particle Analysis 


(1) Inertial Model - (r^ ■ 50y ) 

As discussed in section (B) , inertial effects must be considered for 
sufficiently large particles (> 25vi). On the basis of the conparison used, 
even a 25vi particle had a migration time constant one order of magnitude 
smaller than the diffusion time constant. The question of inclusion of 
turbulent diffusion will be considered in the next section; for the following 
discussion^ diffusion will be ignored. In this section, a particle trajectory 
analysis including particle inertia will be discussed. 

In the Lagrangian coordinates, the momentum equation is 
dv __ _ 

m ^ + 6TTnfVp = eimrU + mg + (4-34) 

where 

U » the mean gas velocity 


mag 


7 " 


2 




U ' 


ji^o m=l 


A normalization is en5)loyed as before, but the base time is an inertial- 
gravitation time. (This normalization is altered in chapter 7 vdiere gravity 
is ignored- -the base time is a characteristic diffusion time). 
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Writing equation (4-34) for the x and y directions in teims of noimalized 
variables gives 


2 r Ax T 

dx ^ in gray — ^ in gray 

dt ^in vis ^ ^^in-vis^ 


(4- 36 (a)) 
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(4- 36(b)) 


In general, one can express (4-36) in four first order differential equations. 
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Equations (4-37) can be numerically integrated if the y integrations (c,d) 
are perfoxmed first. 

A fourth order Runge-Kutta numerical integration system was adopted 
inarching equations (4-37) in time.^^^ The procedure made four approximations 
on the change in v over the time interval At, taking a weighted average of 

“T “ 

these to get v . The change in y over this interval of time is found by a 

“v 

forward Euler approximation. 
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The same procedure is then enployed to obtain and x^ over this time 
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interval with y at midpoint approximated as - - '' 2 • The integration 
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program KDDN4 is listed in Appendix H. The final result is the particle tra- 
jectory* flight profiles for different exit positions from the injection tube. 
Cfigure 4-8). 

The \ise of the flux conservation principle is necessary to obtain density 
infoxmation from these results. In steady state, the spreading of the tia- 
jectory lines is related to the density. Specifically, if the profile tra- 
jectories have an initial separation Ay and nrpact the lower duct wall with 
a separation (ThQ), then 

^ I cn»0^ C4-38) 

y*0 

or after noimalizing 


"o % V ■ " 
o 



Inherent in equation (4-39) is the assumption that the density is uniform 
over the injection tube cross section. 

Inherent in ecpi. (4-39) is the assunption that the problem is two- 
dimensional with no spreading in the transverse z direction. This assumption 
is consistent with using a multi- tube particle injector positioned at several 
z locations across the duct. A transverse slit may be ideal from this view- 
point but would surely perturb the air flow and g^ierate a disruptive wal:e. 
(Figure 4-9 (a)). The z dependence is not critical since all precipitation 
measurements were averaged over the duct width. The effects of turbulent 
diffusion considered in the next section tend to somewhat uniformly spread 
the precipitate after about two duct widths downstream, even for moderately 
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large particles ClOOy) . It was decided therefore to use a single tube 
injector. 

Figure C4*9(^)) depicts a reasonable flux profile i^on exit from the 
injection tube. Quantitatively, the net flux should display the characteristic 
tuiiulent velocity tube profile dependence, i.e.. 


T 

x-0 
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(4-40) 


tidiere 

v„ » the average iniection velocity 
^inj 

r^ ■ the tube radius 

y « the position in, the tube measured from the bottom of the duct 

y^ » the distance to tube centerline from the bottom of the duct 

njj^gjc “ maximum density at the center of the tube 

The net weight of particulate injected during a time t^ in steady state is 

t y *r 

Cinj^ed^ = [ J “ d 

o y«‘r_ 

•^0 o 

and using equation (4-35) 

(tajSted) * Or,} C9 C4-42) 

The flux and weight of precipitate collected on the bottom over an interval 
Ax are respectively: 
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(4-43) 
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Confiining equation (^-44) v/ith (4-42) it follows that the weight collected 
per unit interval Ax becomes 
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Combining the above with equation (4-39) gives 
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(4-46) 


■nnis a loiowledge of the particle trajectories is enough to link the inertial 
theory on an absolute scale to experimental precipitation results. The 
experimental correlation is shown in chapter 6-c. This theory, ignoring 
tuibulent diffusive effects altogether, predicts far too large a particulate 
deposition. 


(2) Hybrid Inertial -Diffusion Model (lOy < < 50p ) 

The results from the time constant analysis indicate migration times are 
an order of magnitude shorter than the duct diffusion time for 55ii particles. 
This was based however on diffusion over -the duct half height. Upon exit 
from the injection tube, the vertical density gradient is quite large near the 
txibe outer radius. Furthermore, chapter 5 reveals that the tpstream diffusivity 


is nearly an order of magnitude larger than the downstream value (.002 
An inclusion of diffusion, effects appears r.ijcfcasary especially in the 
neighborhood of the injection nozzle. The prohleiii then centers on incorpor- 
ating both inertia and diffusion. 

The aim is to accurately predict the density profile and precipitation 
of particles in the diameter range 55-75 microns. The above discussion would 
indicate a pertuibation on the inertia theory analysis is more appropriate 
than a modification of the diffusion small particle theory for this particle 
range. Towards this end^, it appears that the addition of a texm to the 
momentum equation is acceptable because the solution degenerates to the 
diffusion model when inertia effects vanish. A possible solution meeting 
this requirement (bajsically st^erposition of models) is 

dV* n 

* eWpVp - fom-ptr + i€ ♦ - 6imrp ^ (4-47) 

With m » 0, multiplying by density n and dividing by (eTtrirp) gives the flux 
expression (4-1). The diffusion model then follows by taking the divergence. 
Equation (4-47) is in Lagrangian coordinates, and thus the instantaneous 
particle density and gradient must be loiown along its trajectory. The 
evaluation can be made only by incorporating conservation of particle flux. 
This is to be expected since the nature of the diffusion term depends not on 
a single particle’s parameters, but on the proximity of its neighbors. 

One can now proceed along the lines of the previous section. The normal- 
izations mimic those of (4-35) except for the x variables which are normalized 
tc the duct height 2a for ease in evaluating density gradients. 
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As with the diffusion theory of section B, longitudinal diffusion will be 
assumed negligible when ccrnpaxed to convection, l^ith these changes, the two 
noxmalized equations to be solved become 
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The method of integration is identical vdth tlie inertia model. Prior 
to the y integration however, the evaluation of the local density and density 
gradient is formulated for each flux tiibe (from the previous time step results) , 
and this is used as the input for the next time step. The analysis requires 
that, the density be evaluated at every point. As tfie integration is carried 
out in Lagrarigiaii coordinates, it is desired to con 5 >ute the density in. the 
same representation. All particle positions are specified by their initial 
position and time only. Figure (4- 10, a) shows an incremental volume in space 
when t=0, the coiner position designated ? (a,o) = a. By the time t, the 
volume element has changed location and shape as shown in Figure (4-10,b). 

The rate of change of the position vector X with respect to the initial 
volume element edge directions is found by observing the change of these edge 
vectors between the time o and t. 
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Conservation of mass requires that 


(4-50) 


p(a,o) Av (a,o) ■ p(a,t) Av (a,t) 


The initial volume of the differential elements is (A a^ Aa^. A a^) . Ihe final 
volume of the element at time t is just the parallelopiped volume of Fig. 
(4-10,b), i.e. 


4V c?.t) ■ . (|L X II- (4-51) 

This '^oluns is first found by taking the cross product and then carrying out 
the dot product of the volume element edge extentions. This operation on the 
derive of the position vector alone is known as the Jacobies J. Combining 
(4-50) and (4-51) gives 


p(a,t) » p 


(4-52) 


where J = 
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Each element of the Jacobian determinant, e.g. k— , represents the change 

da. 

of position vector C in the kth direction of a particle beginning initially at 
a + Aa . j rather than at a. In a two dimensional evaluation where 

•7 



0 i ^ z 

1 i * z 
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(4-53) 
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then J • 
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The evaluation of the Jacobian and ttius the density at position ' time t, of 
figure (4-11) is obtained from knowledge of the trajectories of three other 
points beginning at position 2,3, and 4 at time t > 0. These elements are 
related to the position vectors as follows: 




^0 



(^31 _ ^41) 

— w: — 
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(4-55) 


where Ax = (v^^ ) At, At being the incremental numerical time step, 
o 

Note that derivatives involving initial changes in vertical position A y^ 
were obtained from the displacement of position vectors to either side of the 
point in question to yield an effective two-sided derivative. This procedure 
was adopted throughout except at the upper and lower trajectories where a 
single vertical derivative was used. 

The final step is to conpute the vertical density gradient from this 
information. The attack employed here is to calculate the true gradient from 
a knowledge of the gradient conponent in two arbitrary directions. Consider 
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again figure (4-11) with the vector rj^i .31 representing the vector from point 

1* tc point 3* and T- representing a unit vector in that direction. 

n'-3’ 

Conponents of the density gradient are found as follows: 


Vn • T 

ri,.3. 
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(4-56) 
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Splitting the gradient into x and y components and performing the dot product 
yields 
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(4-58) 

(4-59) 


Equations (4-58) and (4-59) can be solved to give the vertical gradient 


3 n ^^3 * ^ ^ - ^^ 2 ' ”^1 ’ ^ ^^3 ’ ”^ 1 ' ^ 

W “ (731-/11) (X2t-x^,) - (y2i-yif)(x3,-Xj^,) 


(4-60) 


This analysis is basically good except for the one-sided nature of the vertical 
point difference (Fji^j^,). This one-sided difference was found to precipitate 
oscillations in the computer program. For trajectories other than the upper 
and lower ones, it was better to replace equation (4-56) with one involving 
points 3’ and 4' to approximately represent the gradient at 1', i.e., 


Vn.i 


■4’-3' 


^3* '^4* 
^4»-3' 


(4-61) 


The vertical gradient throughout most of the trajectory region becomes 



(4-62) 
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3n C^ 3 1 “^4 • ^ (^2 ' ^ " C^ 2 ' ^ ^^3* ”^4’^ 

CXji "y4f } C3^2' " ^^2' ^^3’ ”^4*^ 

It was desired to keep the incoming flux profile unchanged from the 
inertial assunption. To also avoid huge gradients at the tube periphery, the 
injection velocity was assumed constant across the injection tube exit area 
and the density given the characteristic turbulent profile, i.e. 



The linking of these results with the deposition on the lower wall proceeds 
exactly as in the previous section (equation (4-46)). The program KDIN7 used 
is listed in Appendix I. 

One final comment concerning particulate spreading by turbulent diffusion 

is in order. Turbulent eddies can never cause movement of particulate faster 

than the turbulent perturbation velocity. It was necessary to put an upper 

limit on at v’ which from cnapter 3 becomes 

n 


V.V 13 

n —NT 


(4-64) 


Furthermore, it was assumed sufficient in the perturbation type analysis to 
keep uniform across the duct core region and setting it to zero one-half 
centimeter from the wall (D^ goes from 0 to maximum value in roughly 1 centi- 
meter) . The theory's insensitivity to the choice of cutoff layer thickness 
shown in chapter 6, justifies this action. 

A comparison of this hybrid theory's precipitation prediction with 
experiment follows in chapter 6-B. The work in chapter 7 concerning particle 
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flight in a jet boundary layer builds on the hybrid diffusion theory. 
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V. DOCUMENTATION OF FIELDS AND FLOWS 

A. Introduction 

The thrust of this thesis Is directed towards predicting the precipita- 
tion of magnetizable particulate In turbulent air flows. Towards this end, 

It Is necessary to quantify turbulent eddy effects and magnetic field ef- 
fects and magnetic field effects. The purpose of this chapter Is three-fold: 

I) to describe the method used for Inducing turbulence In the duct 

II) to document the measurement of the channel turbuelnt diffusion co- 
efficient 

III) to give an account of the permanent magnet array field measurement 

The goal of the first objective was to induce turbulence in the 
duct that would notdecay as quickly as that incited by a screen. It was hoped 
that the design of a velocity profile grid might give rise to a channel in- 
tensity resembling fully developed turbulent flow. The result was in fact a 
turbulent intensity decaying an order of magnitude over the channel length 
(which is considerably better than the turbulence excited by a screen). The 
axial non-uniformity in turbulence was acceptable but necessitated measure- 
ment of turbulent diffusivity at several axial locations. 

The result of the second objective was an experimental measure of the 
turbulent diffusivity's axial and transverse character in the duct. The 
outcome of this testing was the establishment of the credibility of linking 
these measurements to existent fluid turbulence diffusivity literature. 
Specifically, the data supported the decision to use the transverse dif- 
fusivity dependence documented by Davies (chapt. 3 [2]) with the peak mag- 
nitude obtained experimentally along the duct. This course of action was 
motivated by the desire to avoid using transverse experimental results be- 
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cause of the asymmetry they exhibited. The asymmetry was caused by the 
hot wire anemometer probe schroud; the schroud Inhibited fluid flow dif- 
ferently when extended clear across the duct than when extended just past 
the Insertion hole. 

The third objective was to determine the nature of the magnetic field on 
the duct lower surface excited by the permanent magnet structure. The goal 
specifically was to ascertain the relative size of the different harmonics. 

As an aside, It was found that positioning the magnets 1/4" from the duct 
lower surface yielded a nearly sinusoidal field on the duct surface. 

B. Inducement of Turbulence and the Velocity Profile Grid 

This section describes the construction of a velocity profile grid In- 
tended to excite a fully developed turbulent flow. The grid design (an 
Improvement over a plane screen) succeeded In Inducing turbulence decaying 
an order of magnitude down the duct. 

The desirability of a fully developed turbulent flow poses a problem 
in this experiment. Although the Reynolds number for this flow is well 
within the turbulent region (-36,000), the length of duct necessary for 
development of this flow 1s prohibitive. Schlichting [1] states that the 
point of transition at which an instability initiates a growing disturbance 
in the flow has been found empirically to occur when 



-g- * 3.5 x 10° 


(5-1) 


For a 4.3 m/sec flow, this gives a transition length of 2.5 meters. Further- 
more, the boundary layer thickness 5 is given by 
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«(x) - (0.37) X (iiail) -1/5 (5-2) 

Thus the distance x for this boundary layer to grow half the duct width 
(a * 6.35 cm) is 

X ■ “2.6 meters (5-3) 

.37 

It would take over 4.5 meters to Insure fully developed turbulent flow. 

One method of exciting turbulence in the duct is to use a screen or 
wire mesh. According to Baines and Peterson [2] the mesh size of the screen 
and the wire diameter (bar thickness) completely determine the character of 
the turbulence excited. Their results reveal that 5-10 mesh lengths (see 
Fig. 5-1) downstream from the screen are required to insure good flow estab- 
lishment, i.e., homogeneity. It is the screen bar thickness b (wire diameter) 
however, which determines the turbulent intensity decay. The figure shows 
roughly this intensity versus distance downstream dependence in terms of 
number of bar thicknesses. It is desirable to work in the region where 
v'/u = .01, 100 bar thicknesses downstream. This option would however yield 
a very low intensity, and a roughly exponential decaying turbulent level would 
still exist in the expeirmental section. For this reason the option of using 
a nonuniform velocity profile grid rather than a screen was adopted. 

There are at least two requirements that such a profile grid must meet 
to establish a fully developed turbulent flow— 

1) The size of eddies excited by the grid spires must be the same as 
the large scale eddies (L) that would exist in the full developed 
turbulent duct flow. 
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2) The Impedance presented by the grid to the upstream flow must be 
such that the average axial flow velocity U over the duct cross- 
section matches that of the fully developed turbulent flow. 

It Is noted that this discussion does not consider the effect of duct walls 
on the turbulence growth. 

Most work cn profile grids has been done with circular di'cts. A typical 
profile grid for a round pipe Is shown in Fig. (5-2(b)). The grid consists 
of a number of radial spires, the size of which at any radius Is chosen to 
match (I.e., Insure the same velocity of) the fully developed flow. The 
number of spires Is determined 1n meeting requirement (1 ) above. The average 
thickness of a spire (e.g. at half the radius) represents the average eddy 
scale excited; this must be the same order as the turbulent length scale 
for the circular duct flow desired. 

The corresponding grid for the rectangular duct Is shown In Fig. (5-2(a)). 
Each side of the duct has the same number of spires extending toward the cen- 
ter, but each is terminated on the duct diagonals. The same procedure is used 
in determining the width and number of each spire. 

It was concluded after private consultation [9] that the characteristic 
duct Integral ^cale Is roughly a quarter of the duct width. The eddy size 
excited by a spire will be roughly the same as the spire width. This would 
Indicate four spires to a side as being the proper choice, which is what the 
author used. 

The spire width calculation required more work. Frank Durgin [3] has 
studied the excitation of full scale turbulence in pipe flow, and has found 
that the spire width d and the spire spacing i should be related as follows 
(see fig. 5-3): 
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(Ijp)2 . 1 ♦ k -iSZili (5.4, 

where U ■ the average X directed velocity 

Uq ■ the maximum midstream x directed velocity 
« constant 

In circular pipe flow Is about 4, but private discussion [9] has led the 
author to believe k^ * 10 Is more accurate for rectangular ducts. 

To complete the calculation, a formula for flow velocity with position 
is needed. The author has measured the duct flow velocity with position. 

The results for maximum flow speeds of 850 ft. /min. (4.3 m/sec) and 

450 ft. /min. (2/2 m/sec) are shown In Fig. 5-4. The asyirmetry about the center 

line Is attributed to two aspects of the anemometer probe: 

1) At lower y positions, the probe extends across the entire duct 
Impeding flow above the probe tip. At the upper positions leakage 
Is occur Ing through the probe Insertion hole. 

2) The anemometer probe tip Is completely surrounded by a metal 
shroud inhibiting sensitive measurements especially near the 
wall. This is undoubtedly the source of asymmetry exhibited in 
Figures (5-4) and (5-10). 

The author has fitted several curves of the form 

i = (f) (5-5) 

where (n = a constant) to the above data. The best fit appears to occur 
when 0 = 9; these are plotted next to the experimental measurements in 
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fig. (5-4). 

Using the above results* equation (5-4) transforms to 


d 

I 




2 /® - 1 


TiT 


- 1 


1 + 


10 


(5-6) 


Now since A is equal to a quarter of the duct width, d is calculated for every 
position up to the center of the duct. A and a are different for the ver- 

tical sides with rectangular ducts, but (5-6) applies. Thus the author con- 
structed eight identical vertical spires and eight horizontal spires, posi- 
tioning these as in Fig. 5-3 and cutting them along the duct diagonals. 

As will be shown in the next section, the turbulence decay is consistent 
with screen literature predictions if the characteristic screen bar size is 
chosen equal to the thick base (1.25") of the spire. Choosing any other 
thickness of the spire as a representative bar size leads to the conclusion 
that the spires induce a turbulence that does not decay as rapidly as would 
screen turbulence (incidently, 4 screen bars 1.25" thick would block the duct 
completely). The important point is that turbulence is induced and has an 
axial variation. 

C. Turbulent Diffusion Coefficent Determination 

Historically, work in the area of turbulent electrostatic precipitation 
[7] has been geared along one of two limits. One involves the assumption of 
an infinite diffusivity known as the Deutsch model, which implies the fliud 
turbulent eddies dominate all other forces in supplying particulate to duct 
walls where it can be precipitated. The second limit ignores diffusion and 
infers that particles have a large enough inertia that they are unaffected 
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by turbulent eddies. The work by Williams and Jackson [8] assumes a con- 
stant diffusivlty (based on Dhanak's formula —chapter 3) falling Instantly 
to zero a small distance 6 from the walls. The Intent of this section Is to 
build upon the work of Davies (chap. 3, [2]), I.e., to use the transverse 
diffusivlty dependence he documents with the midstream measured duct diffu- 
sivlty. 

The turbulent diffusion coefficient was found in Chapter 3 to be 
Dt = ^ T (5-7) 


— y 

where v ' is 


the time average squared perturbation velocity 


and 


T = 


f v‘lT-t) v*(tT dt 



The long overbar in the Lagrangian time scale T is understood to be an 
average in time t. 

5" 

The turbulent intensity v' was measured with an RMS meter and an ane- 
mometer probe. The velocity probe signal feeds a resistance bridge in the 
anemometer which then in turn supplies a voltage to an RMS voltmeter that 
blocks any dc signal (see Fig. 5-5, a). The RMS voltage is then multiplied 
by a constant to give the turbulent intensity. 

The output from the anemometer wire represents the resistance change 
resulting from the cooling of the hot wire probe. This velocity voltage 
relationship is sometimes referred to as the King relationship and is shown 
in Fig. 5-6. The anemometer signal can be fed into an equalizer to linearize 
the voltage dependence. If the interdependence is linear, then any perturba- 
tion velocity superimposed on the main flow^when filtered through an RMS 
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voltmeter and multiplied by the slope (change In velocity with respect to 
voltage)^ will yield the turbulent Intensity. 

Because an equalizer was not avail ablet a more approximate method had 
to be adopted. Figure (5-6) shows the anemometer voltage-velocity relation- 
ship. As long as the velocity perturbations about the steady mean are small, 
the RMS voltage is related to the average perturbation velocity through the 
local slope. By taking the slope from fig. (5-6) for the average flow at 
every position y In the duct (fig. (5-4)), one gets the proper multiplying 
factor, c, along the duct. The turbulent intensity measured in this manner 
is shown In Fig. 5-7 for mean flow velocities 4.32 m/sec and 2.28 m/sec. 

The measurement of the turbulent time scale Is more involved and re- 
quires a measurement of the energy spectrum, i.e., the turbulent energy in 
Isolated frequency bands. Von Karman [4] and Berman [5], as well as Hinze's 
book ((pp. 165-174) referenced in Chapter 3) use stochastic theory to pro- 
vide a connection between the integral of the velocity auto-correlation in 
eqn. (5-7) and the turbulent energy specturm. Berman, as well as Durgin and 
Fannucci [6], specifically correlate the integral length scale with the fre- 
quency at which the turbulent intensity begins to fall, i.e., 



(5-8) 


where = the maxi urn axial velocity 

f|j = the turbulent intensity break point frequency 

Figure (5-8) shows the turbulent intensity in 20 Hz band widths with 
a main flow of 4.32 m/sec and 2.28 m/sec. The measurements were made in the 
center of the duct in both the x and y directions. The equipment set-up 
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Figure (5-7) Turbulent Intensity Versus Position in Duct 
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Figure 5-8 Turbulent Intensity for 20 hz Frequency Bands 
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1s shown In fig. (5-5(b)). The process is Identical to that above except 
for the addition of the bandpass filter. As fig. (5-8) shows, the 20 Hz 
bandpass test beginning at 10 Hz, revealed, only the characteristic decrease 
in turbulent Intensity with a 5/4 fall off. It was necessary to decrease 
the frequency bands to 4 Hz, centering the first at 4 Hz to observe a con- 
stant intensity and a break frequency. Figure (5-9) shows these results, 
again for measurements In the center of the duct for a 4.32 m/sec flow. The 
axial flow Intensities for each flow exhibit a 19 Hz break frequency at 4.32 
m/sec flow and a 15 Hz break frequency at 2.28m/sec. flow. 

Using equation (5-8) and the work from Chapter 3 in connecting the 
turbulent length and time scales (T = L/U^), the desired turbulent para- 
meters follow; 


4.32 m/sec flow 


2.28 m/sec flow 


f L = 2.7 cm 

T = .00624 secs. 
L = 1.81 cm 

T = .00791 secs. 


(5-9) 


Multiplying these time scales by the respective turbulent intensities of 
fig. (5-7) gives the turbulent diffusion coefficients shown in Figs. (5-10) 
and (5-11). Along with these experimental values are shown the diffusivities 
as predicted by Davies and Dhanak (see chapter 3). Davies's prediction 
based on the duct width Reynolds number gives 
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Figure (5-10) Diffusivlty versus Position for a Main Stream Flow 
of 4.32 m/sec 
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y-cm 

Figure (5-11) Diffusivity versus Position for a Main Stream Flow 
of 2.28 m/sec 



-137- 


» 


.01v(R^ 

®2a 


4.32 m/sec ■ .00147 

(5-10) 

1 2.28 m/sec ■ .00084 


As shown In figs. (5-10,5-11), Davies assumes a rise In diffusivlty to the 
above values In 1/8 of the channel width. Dhanak bases his diffusivlty pre- 
diction on the duct half width Reynolds number as well as the duct friction 
factor. 


D. * .0708 R^ vi/T 


4.32 m/sec * .00205 

( 2.28 m/sec = .0013 


Here it is assumed from Dhanak's empirical results that the friction flow 
factor for a dry air is .011 and .016 for the 4.32 m/sec and 2/28 m/sec 
flows respectively. This value which represents axial dependence only is 
intended to be valid in midstream. The rough agreement of the three is ac- 
tually quite good when one remembers that these expressions are obtained 
empirically for specific flows. 

One further point should be made concerning the generality of this 
approach. The diffu^ivity obtained was based on turbulent intensities 
measured in the flow direction. To apply this diffusivity to diffusion in 
the y direction, one must assume isotropic turbulence, which shall be as- 
sumed henceforth. It was hoped that the turbulent intensities might be 
measured in the transverse directions, but the probe multiplicative con- 
stant is difficult to determine near zero mean flow and the linearity 
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Assurnptlons are questionable. 

The diffusivlty representation adopted In this work uses the transverse 
dependence of Davies with the measured experimental mldstrean: diffusivlty to 
determine longitudinal magnitudes. The hybrid diffusivltles are shown dot- 
ted In figures (5-10) and (5-11). These measurements were taken just before 
the right-most flange of fig. 1-2. One rough order estimate of the dlffu- 
slvlty's axial dependence Is obtained from screen turbulence literature, e.g., 
Baines and Peterson [2]. Given a screen bar size equal to the base spire 
width (1.25”), their results predict a grid turbulence level (~) = .16 and 
a turbulength scale L = .0095 m, 6 inches from the screen. With a 7.62 m/sec 
flow the diffusivlty is 

This would indicate an order of magnitude decay down the duct is likely. 

Fig. (5-12) shows the measured diffusivlty versus axial position. The 
region iiranediately following the grid where the velocity "necks" down ex- 
hibits a calm before a turbulent transition. This data was matched in a 
least squares polynomial fit using the MIT Math library subroutine LSFIT. 

The axial dependence past the 1 ' mark was found to be accurately represented 
as 

* .0546 - .151 X + .202 x^ 

- .0843 x^ - .079 + .0289 x® + .0656 x® 

- .0336 x^ 


(5-12) 
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Figure 5-12 Turbulent Diffusivity versus Axial Position 
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Thls completes the measurement of diffus1v1ty as a function of both x and y. 
D. Magnetic Field Distribution 

As mentioned In chapter 2, the permanent magnet field structure Is not 
sinusoidal, but contains manu harmonics. Given that the field pattern re- 
peats over a length L and begins at zero, a reasonable representation would 
be 


m 




z 

1 


b„ sin 
m 


m' 


2TTmx 


(5-13) 


Multiplying equation (5-13) by sin and summing over 2M-1 discrete X 
steps gives 
(2M-l)Ax, 

(2M-1) Ax,M 

(5-14) 



X*AXi 


(2M-1) Ax,M 

B(x) sin ^ = 21 ZT '’m S’" ^ 


x=Ax-| m=l 


where Ax^ « L/2M 

Using the discrete orthogonality property for sinusoids wherein 


(2M-l)Ax, 

z: 


x=Ax^ 


sin ( 


2Tmix 

L 


) sin ( 


2ttJIx 

L 


•) = 


0 m?«A 

J/2 (2M-1) m=il 


It follows that 


(5-15) 


(2M-l)Ax^ 

I’m = ST 21 Bnet '*> (^> (5-16) 

x=Ax^ 

The normal flux density was measured 1/4” above the permanent magnet 
wave structure for each of three wavelengths— 5.08 cm, 8 cm, and 12 cm using 
a Hall effect prove. For each case, a discrete Fourier analysis according 
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to the above discussion was performed using the computer program FOURIE 
listed In Appendix I. The measured field for the above three wavelengths 
along with the first 9 calculated harmonic components Is shown In figs. 
(5-13) to (5-15). 
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Bj = .0899 = -.000106 

B 2 = -.00156 

B 3 = -.025 

B, = .000267 
4 

B 3 = -.00128 
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VI. Experu»ntal Results and Correlations with Theoretical Models 
A. Small Particle Dynamics 

Figure C^-l) shows the experimental arrangement for the small particle 
tests. The injection nozzle is positioned approximately 3 feet before the 
magnetic structure to insure a uniform initial particulate density over the 
duct cross -section. Several tests appeared to indicate an injection distance 
6" into the settling chamber gave reasonable uniformity. Where possible, 
precipitation was avoided in the local region just after the profile grid 
where flow necking and spurious wakes could cause unwanted side effects. The 
tests were performed for three field wavelengths - 5.08 cm., 8 cm., and 12 
cm. For each experiment, the total weight of particulate injected was re- 
corded, and the amount of material precipitated across the duct width per 
half wavelength axial distance was measured., 

1. Correlation of Measured Parameters from Numerical Models 
It is necessary to correlate the number of grams per half wavelength 
with the density dependence theory calculation. The correlation is based on 
the assumption that particulate at density n^ is injected during a time 
interval t^ seconds at the beginning of the precipitation structure. The 
total weight in grams of particulate injected is 
^0 2a 

CNet Wt.) = [ [ p U n w dy dt (6-1) 

grams ^ ^ 

where 

p = mass density of particles grams/j^^S 

w s duct width-meters 
U = axial flow velocity = 




Evaluation of (6-1) yields 
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(Net wt.) - p(2a)(.9 w n t (6-2) 

grams 


In terms of the vertical particle flux at y ■ 0 , (Ty « q) » the precipitated 
particulate per half wavelength strip in a steady state operation lasting t 
seconds is 


(measured) 

output 


t^ X + X/4 

I j PCiy . w dx dt 

0 X - X/4 


(6-3) 


where x is measured from the beginning of the magnetic precipitation structure. 
The particle flux at y » 0 is in Stokes flow without diffusion so that 



Assuming there exists an average flux over any one axial half - wave strip 
whose magnitude is r^(x)., the measured output (6-3) can be written 


(measured) = p r_ w(i) t 
output 0^0 

where is evaluated at (x + j) 


Using (6-3) this becomes 

( grams of 


\ 


(measured) 

output 




(- ) 


(2a) (.9 W 


o avg. at (x) 


(6-S) 


( 6 - 6 ) 


( 
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2. Correlation of Measured Parameters from Deutsch Model 
The data will be conqjared^to the refined Deutsch precipitation jiiodel 
developed in Appendix B. The Deutsch model analysis correlation proceeds in a 
similar fashion to section 1. The density developed is, 

(6-7) 

and is valid for plug flow U and a single harmonic B field. Including the full 
harmonic field and the actual flow velocity, (607) becomes 




(6-8) 


where is defined in (6-4) 
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Integrating equation (6-3) with ^ q given by (6-4) and substituting the 
total mass injected using equation (6-2) gives 
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3. Observations of Theoretical Deposition Predictions 
The theories of section 4-B are conpared by examining the percentage of 
injected particulate collected per half wavelength on the lower duct surface. 
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Figure (6-2) shows the deposition predictions (with 4^ u particles. B cm 
field wavelength) for the four light particle theories— Deutsch, non causal, 
causal perfunctory, and causal fundamental as presented in chapter 4. The 
numerical theories are computed with the spacer grid representation discussed 
in section 4-C using 25 transverse points and 80 axial steps. Studies re- 
vealed the accuracy of numerical solution with such a grid assignment is 
questionable for two reasons, figure (6-2) is useful for comparing theories 
only in a cursory manner. 

The first problem to be considered is that of numerical convergence. 

The causal theory downstream (7.5 cm) deposition prediction decreases 505S 
as the grid size is increased from 15 to 40 transverse points (axial » 80). 
Considerations of grid stability indicate the axial spacing Ax should be less 
than or equal to (D|.U). This criteria sets the number of axial grid points 
into the thousands. The numerical density solution was calculated in the 
limit of constant diffusive and convective terms only, and compared 
to the analytical solution. When using 95 transverse points and 1000 axial 
steps the numerical solution demonstrated reasonable convergence and dif- 
fered from the analytical solution o ily in the third significant figure 
(<1% difference). 

The second problem of numerical accuracy is more subtle. The diffu- 
sivity and magnetic force in the lower portion of duct near the laminar 
sublayer change very rapidly. These facts added to the uncertainty as to 
the natureof the sublayer physics make it clear that the lower portion of 
the grid solution at the wall is the area of greatest uncertainty. The 
solution outlined in section (4-C) couples the migration fluxes to the 
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density at the lower grid point (for each axial location) to arrive at de- 
position on the wall, although the net profile is relatively unaltered. 
Considerable changes in predicted densities occur for the lower transverse 
grid point locations when the number of transverse points is increased 
further. 

An alternative to predicting deposition rates which effectively in- 
tegrates through the region of numerical diffuculty, Involves applying thi 
principles of mass conservation to the numerically calculated density pro- 
files (for which the analytical model gives credibility). Specifically the 
transverse flux deposition on a strip Ax equals the difference of upstream 
and downstream axial fluxes. The procedure of using the density profiles 
of section (4-C) in this manner to arrive at deposition predictions is 
described in Appendix K, and the corresponding refined causal fundamental 
FORTRAN program KDFFl is listed in Appendix L. 

Figure (6-3) shows the causal and refined causal deposition predic- 
tions (grid 95 x 1000) along with the Oeutsch theory prediction for 8y 
particles and an 8 cm field wavelength. The increased grid accuracy re- 
veals the error of the unrefined model. Increasing the grid accuracy from 
(95 X 1000) to (500 X 2000) gave only a 5% decrease in deposition pre- 
diction (the 95 X 1000 grid assignment will be used for the remainder of 
this chapter). The refined perfunctory prediction exceeded only slightly 
the fundamental causal preductions (<1% difference). 

The Oeutsch model, with an infinite diffusivity always predicts greater 
precipitation (except at the initial axial locations) than the finite 
diffusivity models and exhibits a less rapid initial decay. This observa- 
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tlon Is linked to a property common to all the theories concerning the 
shape of the deposition prediction. Mass conservation dictates that the 
predictions must have the same integral (i.e.» net area). Thus, steep 
initial deposition predictions must be balanced by slow downstream decays. 
Theoretical predictions are limited in how they can differ; the principle 
difference occurs in the nature of the initial deposition decay. 

As pointed out in fig. (4-6), the diffusivity undergoes a rapid change 
in slope at the laminar sublayer. An attempt to simulate this change in the 
causal fundamental model (in which 0^ goes to zero Just above the sublayer) 
was made by increasing the magnitude of the slope by a factor of 10 for the 
last transverse core grid point only. This assignment is in fact an arti- 
fice, but may represent the diffusivity dependence more realistically. The 
numerical deposition prediction was not sensitive to this alteration; leav- 
ing the diffusivity slope unaltered results in a 0.5% lower deposition pre- 
diction (at X = 7.5 cm). The factor of 10 increase in diffusivity slope 
will be inserted into the last grid point assignment for the remainder of 
this section. 

The refined causal fundamental model displayed a much smaller sensi- 
tivity to particle size than did the Deutsch theory (4^ y tests were 13 to 
20% lower at upstream locations, negligible downstream). It is expected 
that small downstream precipitation predictions should exist. Magnetic 

n 

migration is extremely dependent on particle radius (-r')i at upstream 
locations, magnetic forces remove particulate in the vicinity of the wall at 
different rates dependent on their size. After particulate is removed 
from this magnetic interaction region, further removal is more dependent 



on turbulent diffusive effects which are independent of size. The only 
parameters Important in the Deutsch model spring from migration forces which 
depend on the particle diameter squared. 

Discussion of the representation of the field harmonics and their 
effect on theory predictions is in order. Numerical predictions differed 
little (approximately U for X * 12 cm) when only one, rather than nine 
harmonics (chapter 5) were used to represent the field. The nine harmonic 
tests cannot be treated by superpositions, i.e. , more harmonics yield a 
higher flux toward the wall, and thus greater precipitation. The actual 
field being more uniform in sections has a smaller field gradient, and thus 
a smaller force. The fundamental harmonic field computer test predicts 
greater precipitation. A nine harmonic field representation is used through- 
out this chapter. 

Figure (6-4) shows the effect of removing the magnetic field altogether. 
The deposition decays slowly downstream. Such a drop in precipitation is 
not surprising since the magnetic force is 81 times larger at the wall 
(X = 8 cm) than the gravitational force. Only a dusting of the duct lower 
wall was observed experimentally when all magnets were removed. Reentrain- 
ment made a true measurement of gravitationally deposited particulate dif- 
ficult. 

Figure (6-5) shows the theory sensitivity to the magnitude of D^; 
turbulence at larger diffusivities keeps particulate dispersed evenly over 
the upper portion of the duct and sweeps more particulate to the wall where 
it can be precipitated. The numerical solution should converge to the 
Deutsch prediction as diffusivity increases. The Deutsch precipitation 
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Figure (6-4) Refined Causal Fundamental Precipitation Results; Magnetic Versus 
Gravitational Migration, X = 8 cm, 4 1/2 p Particles 
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predictlon exceeds the causal model prediction further downstream, because 
particulate Is not swept up to the top of the duct by eddies. In the limit 
of zero diffusivity, particles decay along trajectories; in section (2-C) 

It was shown that particles In trajectory precipitation are removed more 
quickly than the Deutsch model predicts. The point to emphasize Is that 
precipitation at upstream positions Is quite sensitive to the value of dif- 
fusivity. 

4. Numerical Density Profiles 

Figures (6-6) to (6-8) show the density profile transverse dependence 
at four axial locations for 5.08, 8, and 12 cm field wavelenghts respectively 
with 8y particles as predicted by the fundamental causal model (95 x 1000 
pts). The density rises from zero at the top of the duct, peaks below the 
center, and decreases toward the lower wall. The greater the penetration 
of magnetic forces Into the volume (i.e. at larger wavelengths), the greater 
Is the decrease In density profile magnitude downstream due to particulate 
removal. This is confirmed in figs. (6-6) to (6-8). 

The lowest density plotted is actually that calculated at the laminar 
sublayer (y/2a = .002). The final wall density ranges from .92 to .95 of 
the density at the sublayer (at y = .002*2a cm) as predicated by the trajec- 
tory model of Appendix B. 

The density gradient Is greater in the sublayer than in the core (see 
fig. (4-2)). The perfunctory model Imposes a zero gradient at y = A, and 
as explained in chapter 4, it does so by assuming the density just below the 
core region to be the same as that above. This nunierical assignment is 
certainly ad hoc, since the density does not increase in the sublayer. 
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Figure (6-6) Refined Causal Fundamental Density Profiles with Dlffusivlty and Magnetic Force Decay 
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Figure (6-7) Refined Fundamental Causal Theory Profiles with Corresponding Dlffuslvity 
and Magnetic Force Decay, X = 8 cm y Particles 
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Figure (6-8) Refined Causal Fundamental Density 
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The splicing on of the trajectory model solution distorts what the core 
diffusion analysis assumed should exist below y ■ A. The diffusivity of 
fig, (4-6) undergoes a power law decay to zero at the boundary layer, not a 
discontinuous decay. This implies from Appendix D that the sublayer density 
gradient be much larger than the core gradient. The fundamental model pre- 
dicts this condition without imposing a difficult to justify boundary con- 
dition, and inherently represents a more realistic diffusivity by keeping 
D|. finite in the core region. 

The fundamental model deserves the title "fundamental" for another 
reason. As the degree of turbulence diminishes, it is meaningful to think 
in terms of trajectories. In this limit the fundamental model imposes 
boundary conditions only where particle trajectories enter the region of 
interest. 

Included with the plots of (6-6) to (6-8) are the diminishing 
diffusivities at the given axial locations as well as the exponential mag- 
netic force decay. The results indicate magnetic migration is relevant 
over only 1/4 of the duct. The density profile resulting when magnetic 
forces are eliminated entirely is shown in fig. (6-9). The results con- 
firm the proposition that the densities are not altered significantly over 
the upper three fourths of the duct. The smaller migration over the lower 
1/4 of the duct yields a smaller deposition despite the large wall density. 
The densities in the uniform migration (no D^) sublayer do not decay. 
Finally, it is interesting to consider the density profiles in the limit 
of no turbulence, i.e. = 0. The density (fig. (6-10)) can only decay 
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Figure (6-10) Transverse Density Profiles with Corresponding Magnetic Force Decay, No Diffusivity 
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over the lower portion of the duct where magnetic forces are significant. 

(No density decay occurs over the region of divergence-free force where 
gravity dominates.) The core region of the duct should have uniform density 
dropping discontinuously to zero near the upper region. The discontinuous 
density jump to zero represents the point abjve which all trajectories can 
be traced back to an entry point on the upper wall where the density is 
zero. All trajectories below the transition point have their origin at the 
entry plane x*0. 

An examination of the experimental results and their correlation with 
theory will follow a discussion of the particle size distribution. The non- 
causal model will not be considered in this comparison; its predictions are 
for rcaciily seen reasons, too low. The perfunctory model's predictions 
differed only slightly from the fundamental model. For this reason and the 
above discussion concerning difficulties in justifying the zero gradient 
boundary condition, the perfunctory model will not be compared with date 
either. 

5. Particle Size 

The iron powder used in the diffusion tests ranged in diameter from 1 
to 26 micrometers. It was necessary to know to what extent the particles 
acted as agglomerates in flight. A microscope analysis (histogram) of the 
particle agglomerate size spread is shown in fig. (6-11). (The tests was 
prepared by blowing powder through the injection tube and allowing it to 
settle on an oil coated microscope slide.) The scope base particle size 
was about 4 micrometers. With the intent of determining more of the 
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particle's size behavior in flight, an Anderson Impactor Test (which balances 
a particle's inertia forces with drag forces) was conducted. The results 
(fig. (6-12)) indicated that then equivalent aerodynamic diameter represent- 
ing particles dispersed essentially as in the magnetic collection experi- 
ments ranged primarily from 4^ to 8 microns in diameter. The microscope 
test procedure would favor collection of heavier particles; this may be 
the reason for the disparity. 

6. Theoretical Predictions Compared to Experimental Deposition Data 

The Deutsch and fundamental causal theory predictions are compared to 
experimental collection results at field wavelengths 5.08, 8 and 12 cm 
tests in figs. (6-13) to (6-15) respectively for particles as described in 
the previous sections. Included with Fig. (6-15) is the fundamental causal's 
prediction of deposition of Zu particles. The data supports the causal 
model 's steeper decay over the Deutsch model's gradual decline. 

The Deutsch model prediction increasingly departs from the data with 
decreasing wavelengths. The micron theory results shown in fig. (6-7) 
reveal that th^ Deutsch model still predicts a high particulate deposition at 
a 5.08 cm field wavelength, and a slow axial deposition decay. It was noted 
in chapter 2 that the Deutsch model is based on two major assumptions— 

(1) Fluid turbulence is infinitely capable of supplying particles to 
the wall . 

(2) The migration forces (primarily magnetic) hold the particle to the 
wall; their interaction in the duct is small. 

Assumption (1) is always optimistic. At short wavelengths where assumption 
(2) is accurate, the Deutsch prediction should be high because of assumption 
(1). At the larger wavelengths (with the magnets spaced further apart), 
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Figure (6-13) Refined Fundamental Causal Model Predictions with Deutsch and Experimental 
Data, X = 5.08 cm. 
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Figure (6-14)Fimdamental Causal Refined Model Predictions with Deutsch and 

Data, A = 8 cm. 
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Figure (6-15) Refined Fundamental Causal Model with Deutsch Predictions and Experimental 


Collection, X = 12 cm. 
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forces are larger at further distances from the wall, and thus augment the 
precipitation process (the optimum wave number Is approximately the recipro- 
cal duct height). The error of both assumptions then emerges In a trade off 
at larger wavelengths where the Oeutsch model departs less from data. 

7. Discussion of Results 

Large initial deposition is expected since magnetic forces remove parti- 
cles immediately at upstream locations In the magnetic field dominant region. 
Larger upstream theory depositions are in part caused by the step Input in 
Initial density (high density edge gradients; sine wave input yielded lower 
Immediate initial deposition). Dividing the most uncertain parameter D^ by 
10 (which makes upstream D^ = 2*(literature value) -developed flow) improves 
agreement with data. Altering the second uncertain parameter r (fig. 6-15) shows 
significant alteration in initial deposition also. In the light of these 
comparisons, model credibility along with a useful degree of accuracy in 
predicting precipitation and density profiles exist. 

Three dimensional effects are the most significant effects deleted from 
the theory. The true velocity representation as well as magnetic field ini- 
tial and edge gradient effects should enhance precipitation. Fruitful 
future research may be directed toward understanding the effect of precipi- 
tated dendrite augmented tip effects. Certainly one starting point 
for determining the magnitude of this effect is by determining whether the 
rate of deposition is altered by the loading of already existant precipitant. 

This effect of augmented gradient effects is a concerted interest to re- 
searchers in mineral beneficiation and particulate removal. 



Of the effects discussed above* those related to wavelength are 
the enhanced gradient effects. Specifically the augmented field gradients 
due to magnet edge effects should become more prominent at large wave- 
lengths. Furthermore, non-real i Stic representation of the transverse 
dependence of axial velocity may be wavelength dependent as well. Too 
large a velocity was adopted throughout the core analysis except at the 
center of the duct. Slower moving particulate would in reality have more 
interaction time with the field forces. This effect would be more sig- 
nificant when the force penetrates into the volume (i.e. at larger wave- 
lengths). 

It is hypothesized that those three-dimensional effects are not 
as important as accurate representation of the axial and transverse dif- 
fusivity dependence. Despite the lack of these refinements, the refined 
causal fundamental model has a useful degree of accuracy in predicting 
sedimentation in particle flows where inertia is unimportant. 
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B. Large Particle C50-100u) Analysis - Inertial Hybrid Model 

Figure (6-16) shows the physical arrangement for the large particle tests. 
Lhlike figure (6-1), the injection tube extends well into the duct, usually 
14'* past the profile grid. The number of grains injected as well as the amount 
precipitated per half wavelength is recorded for each test. Except for an 
isolated case, the injection velocity was adjusted to match the main duct flow 
in order to avoid secondary fluid mechanical effects due to non- isokinetic jet 
mixing. The jet was found to cause the diffijsivity to increase only slightly 
in the duct. The new diffusivity profile was fitted with a least squares 
polynomial fit and used in the hybrid theoretical model discussed in Chapter 
4, section D. The first nine harmonics of the magnetic field structure \ised 
in each test were used for the magnetic force calculation. The position of 
the injection tube varied roughly from one -third to two -thirds of the total 
duct height. 

(1) Diffusion- free Model Results 

The first test (figure 6-17) reveals the inaccuracy of the trajectory 
precipitation model (Appendix B) and the inertial model (chapter 4, section 
D) for 60 micron particles. Plotted are the particle paths from the top and 
the bottom of the injection tube for the two trajectory models along with 
the percentage particulate collected per half wavelength (6 cm) . Particulate 
was injected isokinetically at 1500 ft/min. The inertia- free model predicts 
precipitation far too early primarily because of the immediate vertical 
viscous -limited velocity upon injection. The particle's inertia keeps the 
particle in flight longer and predicts precipitation in nearly the same region 
as the observed experimental collection peak* However, the inclusion of 
diffusion is clearly necessary if the trajectory spread is to have any 
correlation with the experimental observations. 
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Figure (6-17) 60 Micron, 12 cm Wavelength, 1500 ft. /min. Injection Trajectory Models with Experimental 
Collection; the left axis pertains to trajectories, while the right pertains to data deposition 
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Figure C6-1S) shows the same test with 100>180 micron particles. The 
experimental collection spread is not as large but again exceeds the inertial 
model prediction. The peak collection again occurs near the inertia model's 
iirpact prediction. The viscous drag dominated model also predicts precipitation 
much too early. The percentage collection consistent with the above trajectory 
models was not shown because of its large value , being at least an order of 
magnitude larger than the observed collection for the 60u case. 

(2) Hybrid Inertial -Diffusion Results 

Figure (6-19) compares the e:)q)erimental precipitation with the hybrid 
inertia-diffusion model prediction. As with the diffusion tests, the per- 
centage particulate (of the total injected) collected per half wavelength 
is the precipitation ordinate. For this 12 cm wavelength, 70 micron particle 
test, the hybrid theory is seen to predict a conparable peak collection 
shifted downstream somewhat from the observed peak collection. The particulate 
spread is much larger now but still smaller than the experimental prediction. 

For illustration, the inertial theory prediction is shown at the top of the 
graph quite compressed in spread with a predicted precipitation an order of 
magnitude larger than the test results. 

The results of 5.08, 8, and 12 cm wavelength tests for 60-70 micron 
particles are shown in figures (6-20) to (6-22). In all three cases the 
inertial diffusion hybrid model predicts larger precipitation with a peak some- 
what downstream of the observed peak. Experimentally only 1/3 to 1/2 of the 
material injected was collected over the approximate four foot working section. 
The theory appears to always predict a short total collection distance, but 
this could only be thoroughly investigated with a longer duct working section. 

(3) Particle Size Considerations 

Particle size sensitivity is pertinent in all models. The iron particles 
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Figure (6- Ig) 100-100 Micron, 12 cm Wavelength, 1300 ft/mln. Injection Trajectory Models with Experimental 
Collection; the left axis pertains to trajectories while the right pertains to data deposition 
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Figure (6- ‘20) 5.08 cm Wavelength, 1000 ft. /min. Injection, 70 Micron- BybriA Model vs. Data 

51% Particulate collected In experiment 
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used in these tests were sized by sifting through several progressively 
smaller screen meshes. The smaller size particles used above were caught 
between 38 and 45 micron mesh screens, the hypotenuse of this size square 
opening being 54 and 64 microns respectively. The largest spherical particle 
must indeed be between 38 and 45 micron; most particles were however anything 
but spherical. Thus, ellipsoidal and elongated particles much heavier than a 
45 micron sphere could surely exist in a batch. Observation of the particle 
group under the microscope revealed 60 micron to be the representative 
diameter of the 38-45 micron batch. 

The particle size sensitivity of the hybrid model’s particulate deposition 
prediction for the 8 cm wavelength case is shown in figure (6-23), As with 
the diffusion case, the size dependent parameters are proportional to the 
particle radius squared, a 20 micron spread appears to be quite significant 
indeed. The interesting result here is that the smaller particles with less 
inertia have a smaller deposition per half wavelength (over most of the 
precipitation region) and a peak deposition i^stream of the larger ones. 

Figure (6-24) shows the hybrid model's deposition prediction for 40 micron 
spheres in this 38-45 micron sifted batch for a 12 cm wavelength field. 

Figures (6-25) and (6-26) show the equivalent conparisons for a 5,08 cm test 
45-54 micron batch and an 8 cm, 38-45 micron batch. 'Ihe upper trajectories 
in the 8 cm test are subject to question because they involved collisions with 
the upper wall. 

The improved agreement leads one to believe that larger oblong particles 
behave in flight more like smaller spherical particles. This idea would 
certainly have credibility in terms of viscous drag forces. The streamlining 
of the heavier oblong particles and resulting lower drag would play against 
the larger gravitational and magnetic forces. The microscope, two-dimensional 






|Ti 


Figure (6-24) 38-^5 Micron Batch9l2 cm Wavelength* 
40 Micron Spherical Prediction 
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Figure (6- 25) 38-45 micron Batch, Comparison of 60 and 40 mlpron Theory Predict Ions, X=8cm 
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slide picture is unfortunately lacking in its ability to reveal the total 
particle shape. 

(4) Boundary Layer Thichiess Sensitivity 

For conpleteness f the question of model boundary layer thickness sensitivity 
must be addressed. The deposition predictions of the hybrid model at an 8 cm 
field wavelength with the boundary layer thickness varying from 1/4 to 3/4 cm 
are shown in figure (6-27), Unlike the diffusion case study, the results are 
relatively insensitive to thickness choice, varying only slightly at downstream 
duct positions. This is understandable because of the much greater role the 
migration forces should play with larger particles, especially the magnetic 
force (which grows exponentially towards the wall) . 

(5) Non- Isokinetic Injection 

In all the precipitation tests, an effort was made to inject particles 
as close to the ambient duct flow velocity as possible. Figure (6-28) shows 
the collection results and theory predictions for a non- isokinetic test where 
the jet injection (2300 ft. /min.) was more than twice the duct velocity 
(925 ft./min.). The working duct section was not long enough to observe the 
con^jlete deposition pattern; less than half of the injected particulate was 
collected. The theory deposition prediction was based on a 70 micron particle 
representing a 45-53 micron grid batch (which, based on the above comments, 
is subject to question). Nevertheless, the results indicate a substantial 
difference between theory and experiment. 

As figure (6-29) points out, there exists a region at the top and the 
bottom of the jet in which momentum is being exchanged between the two fluid 
streams. The thickness of the layer grows approximately as 
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Figure(6-2 7)Hybrid Model Boundary Layer Thickness Sensitivity 
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Figure ( 6 - 28 ) Non-isokinetlc Particle Injection* 8 cm Wavelength* 2500 ft/min, .injection, A5-53 micron batch 
Hybrid Theory and data; 39% collected, 70 micron 




Figure (6-29) Growth of Turbulent Interchange Region 
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The fluid within this region 6 is highly turbulent. Measurement of the 
fluid turbulence from such an injection showed little significant increase in 
diffusivity, but that may be attributed to the insensitivity of the anemometer 
probe and equipment. 

(6) Field Harmonics 

Finally » the inclusion of additional field haimonics has negligible effect 
on the precipitation level. Light particles are supplied by a relatively, 
strong diffusion force to the wall. The heavier particles depend on the 
magnetic force and gravity to get to the lower wall. As the fluid field 
decays exponentially into the volume (the harmonics to a greater extent 
their omission or addition changes little. 

(7) Conclusion 

The effect of adding diffusion to the momentum equation is to spread the 
particles out. The technique of marching in time in a Lagrangian frame, 
calculating density locally at every point in space is applicable to any 
geometry provided knowledge of the turbulent diffusivity is available. The 
diffusive term becomes more important when steep density gradients exist 
Ce.g. at injection tubes). The next chapter will use the hybrid diffusion 
to analyze particle flight over a flat plate. 
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VII. Particle Flight Over Jet Fuselage 
A. Introduction 

The possibili"*^ of ducting particles with a pexmanent magnet structure 
will be considered in this section. The objective is to retain a sibstantial 
part of the particulate within the boundary layer as it evolves. A first 
atteupt at this is to apply the hybrid inertial precipitation work of Chapter 
4 to the problem. 

The hybrid diffusion heavy particle computer model substantiated in 
chapter 6 will be the tool for analysis of this theoretical study. The flow and 
boundary layer development over a flat plate will be used to represent to 
first order the aerodynamics over an aircraft fuselage. All airstreara com- 
pressibility effects will be ignored. 

B. Fluid .Mechanics and Diffusivity 

Any analysis requires first, knowledge of the boitidary layer thickness 
and diffusivity. The experimental evidence (Schlichting-see reference Ch. 5) 
indicates that if the boundary layer becomes turbulent at axial position 
x=Xq, its thickness will increase as 


5 = 



1 

T 


C7-1) 


where v = local atmosphere kinematic viscosity and 

Uq * anbient air stream speed relative to the jet 
The transition to turbulence occurs when 

r U X . 

3 X 10^ < -~ < 3 X 10° 


(7-2) 


For a passenger jet cruising at 40,000 feet, air speed 500 MPH, the onset of 
turbulence occurs about 24 inches after the leading edge and the boundary layer 
thickness is .68 meters, 10 meters downstream. 


The tuxbiilent diffusivity depends not only on the axial displacement 
along the fuselage, but also on. the noimal position in the boundary layer. 

Hinze (see chapter 2 reference, p. 645) shoi^s Klebanoff's data for the eddy 
viscosity’s distribution acrc^c the boundary layer. Equating the eddy viscosity 
with the tuxbulent diffusivity (as discussed in chapter 2) yields the diffusivity 
profile shown in figure (7-1) . A least squares polynomial fit was made to 
represent D^. The results are 

- .037Uq6 [-.00103 + .291 (^) + .007 (^) ^ 

-1.36 (^) ^ + .418 ^ +2.123 (J) ^ (7-3) 

-.277 ( 1 ) ® ■ 1.5S6 (^) 8 . 1.20 [}) 

Typically, injection would occur a few millimeters vertically into the 
boundary layer. As shown in figure (7-1) , the diffusivity near the surface 
becomes negligibly small at downstream positions where 6 is large. The 
diffusivity can be approximately represented as 

= (.037Uq 6)(.07) sin (7-4) 

Using tl\e expression for 6 from (7-1) , the axial position for maximum is 
found by setting the derivative of (7-4) equal to zero. Tiie lesult is 

tan ^ (7-5) 

The graphical solution of (7-5), shown in figure (7-2), reveals that the only 
meaningful solution occurs where ^ goes to zero. The actxial vertical 
dependence (7-3) of D.^. goes to zero faster than sin|^ indicating that for a 
given vertical location, diffusion will have a more pronounced effect at up- 
stream (smaller boundary layer thickness) locations. 
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C. Typical Injected Density Profiles 

Figure (7-3) shows coii 5 )uter calculations of particle trajectories for 

particulate injected isokinetically through an Siran diameter tube into an 

aerodynamic boundary layer 6ram away from the plane skin. Injection for this 

calculation was set at just over 1 m downstream of the leading edge where the 

boundary layer thickness had grown to 7.3 cm. The calculation was identical 

to the hybrid inertial diffusion model developed in chapter 4 except that 

gravity was ignored here. The calculation was perfomed assuming use of 40 

micron spheres, an eight cm wavelength field structure, and an ambient air 

*?p^ed of 500 MPH. The strength of the magnetic field was chosen to be conpa- 

rable with commerically obtainable permanent magnets --the field strength on 

the aircraft skin 1/4" from the magnet being about .095 web/m , Iron powder 

was chosen as the injected particulate, but any large u material of density 
3 3 

~7 X 10 kg/m would give equivalent results. A listing of the program KD747G 
used is listed in Appendix J. 

The effect of the small layer thickness and thus large initial diffusivity, 
is seen to force the upper trajectories away from the plane structure. The 
lower trajectories are, however, precipitated rather quickly. It is not until 
the lower trajectories begin to precipitate and thus lower the densities 
among the upper trajectories, that they begin to drop as well. This i^ard 
diffusion makes it quite difficult to collect the last 5% of particulate. In 
the above calculation, the upper trajectory was only beginning to curtail its 
upward ascent 47 meters downstream, 7.5 cm above the skin. Because of the 
diminished effects of diffusion downstream at the same injection height, 
precipitation of all the particulate is easily accon 5 )lished (e.g. with ^injection 
= 10m, all particulate is collected in 3 meters). In any event, with an initial 
density injection profile characterized by a (y/a) ' dependence, ducting is 




Figure (7-3) Particle Trajectories Over Jet Fuselage 
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not possible; either the particulate will be precipitated or lost. 

D. Possibility o£ Quasi-Stationary Profiles 

The question that mist be considered is whether or not there exists an 
injection density profile that would yield a steady or quasi-steady positioning 
of particulate in the boundary layer. This would require a balance between 
the Inward diffusive flux and the downward magnetic migi*ation. From the 
results of chapter (4-D) , this would require 


6iTTirD, 
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t ^ 3n 
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‘rrr* 




mB^ 0 

1 m 


(7-6) 


or with a fundamental harmonic field only 




(7-7) 


Equation (7-7) could be integrated if "he density at some location y 
were loiown. Although the assignment of n=0 at y*6 is reasonable, also goes 
to zero there and the evaluation of the first term above becomes a problem. 
Again approximating as (.037 11^6) (.07) sin|^j. The analysis proceeds by 
assuming 
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Equation (7-7) gives 


ake'2’^ = C6w)(.037 U^S) (.0?) [| (S-y)j ^ J C7-9) 


or 


ake 
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The constant in front of the exponential is equal to about 82 for parameters 
characteristic of the above text. The e^qionential power 2k6, however, is 41, 
indicating that p is a vexy small number. The density is nearly uniform, but 
drops sharply at y~6 in the neighborhood of the boundary layer exterior. Desig- 
nating the density at y»6_ as n^, equation (7-7) can now be integrated to give 
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idn 
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ake 
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dy 


5(6TTnr)(.037U^).07 sin ^ 


(7-11) 


Reversing the order of integration, (7-11) becomes 


- (y- • 


ake 


•^^4) d(^) 


(6imr)(.037UQ).07 sin ^ 


(7-12) 


Equation (7-12) was integrated numerically beginning 10 ft. downstream 
from the transition to turbulent instability where 6=.26m for an 8 cm wave- 
length field. The calculation was quasistationary in that it assumed a constant 
boundary lay" ■ thickness. The resulting density profile is shown in figure 
(7-4) alow ith the necessary sharp cutoff dependence near y=6. 

The .rge density at y=0 indicates a problem occurring in this region. 
There is no way to terminate the density in this region except as a step at 

y=0^ from n*0 to a very large n almost instantaneously at y=0+e. At this 
point, the search for a possible steady profile becomes merely a mathematical 
exercise, since the profile in this region implies a very rapid precipitation 
ratio. 


It is interesting to compare the spreading of two density profiles under 
the same conditions. The trajectory profiles with isokinetic injection of 
particulate 3 meters after the transition to turbulence, under the same condi- 
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tions as the trajectories in figure (7-3), are shovm in fibres (7-S) and (7-6). 
TTie characteristic turbulent profile (y/a)^^^ dependence is assumed in fig. 
(7-5), while an initial density profile similar to that in figure (7-3) is 
assumed in fig. (7-6). The calculated densities are plotted for various 
locations downstream. The results indicate the skewed profile case does 
indeed keep the mid to upper trajectories more horizontal for a longer axial 
distance. The precipitation of the lower trajectories occurs within the same 
distance as the noimal density profile case however. Downward diffusion 
dominates at these lower trajectories since no particulate exists between the 
plane skin and the lower part of the injection tube initially. The density 
plots reflect just how quickly the effect of downward diffusion and migration 
of the lower trajectories propagates upward. The author's conclusion from the 
above analysis is that ducting of magnetizable particles in a turbulent boundary 
layer with a permanent magnet structure is unfeasible. Ninety percent of the 
particulate injected 6 mm above a flat surface in a 500 MPH airstream can be 
contained end collected within 4 m. The last 10% if collected, will be spread 
over the following 4 m stretch. 
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VIII . Concluding Remarks 

A. Particle Entrainment over Jet Fuselage 

The process of entraining magnetizable particulate in a boundary layer 
may now be examined in the light of the chapter 7 results. If particulate is 
moving in near synchronism with the jet air speed, and collection of particulate 

f 

and reintrainment at the front of the craft is in effect, it is reasonable to 

assume particulate undergoes one full cycle when the jet has traveled twice its 

length mth respect to the air at ’‘infinity". Th^ density tests of chapter 7 

indicate that at best 4% of the injected particulate will be lost per cycle. 

(This assumes that particles collected prematurely are shuffled along the wall 

while still interacting with the boundary layer in some positive drag reducing 

manner.) Thus with these figures, a 100 ft. craft undergoing a 10 mile flight 

**■ C10)C5280 

(and hence round trips) 100 [1 - (1 - .04) ] » 99.9979 per cent 

of the original particulate would be lost. 

Assuming the craft has a surface area of 4700 sq. ft. and specifying that 

1 gram of powder be exposed to 2 sq. ft. of the craft skin at any instant, it 

would be necessary to begin the flight with 1.13 x 10^ kg of powder stored on 

board! Clearly this is an unacceptable demand and an alternative in idiich 

constant interaction with particulate over the skin must be sought. 

B. Travelling Wave Interaction with Particulate 

A travelling wave structure mentioned in the introduction has been 
examined qualitatively as a means of shuffling particulate along the skin in 
an aerodynamic boundary layer. The experiment shown in fig. (1-2) involved 
injecting particulate just under the duct iqjper edge; particles were captured 
by the field, held against gravity, and shuttled downstream opposite to the 
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direction of the travelling magnetic field. 

Figure (8-1) shows an inverted picture of the field region just above the 
wave structure with a right traveling wave, an observer at point P encounters a 
counter clockwise rotating field. Non- spherical agglomerates experience a counter 
clockwise torque, in such a field resulting in a reverse walking motion on a sur- 
face adjacent to the structure. 

Figure (S-^i depicts a typical instantaneous picture of agglomerates on the 

duct i:^>per wall. Particles generally fom agglomerates in the magnetic field. 

% 

Most of the agglomerates will cling to the end of the field structure. Agglom- 
erates over the body of the structure actually bounce along the surface toward 
the downstream field edge. Particles and agglomerates gravitate to some degree 
toward the side edges as they propagate down the duct because of the edge effect 
gradient. 

The flight of a single agglomerate is shown in fig. (8-3) . The agglomerate 
literally walks end upon end; at low frequencies- (1 IIz) this effect becomes clear 
as long agglomerates (1/4") walk along the duct. 

Agglomerate speed is a function of frequency and agglomerate length. Two 
steps are taken every cycle, each equal to the agglomerate length. Agglomerate 
size has a strong dependency on wave structure frequency, there being a gradual 
decay in agglomerate length between 1 Hz and 40 Hz. Above 40 Hz particle 
agglomerates are small ( 400v) . Inertial effects become even more important 
above 60 Hz where motion is impeded. Particles are observed to form into very 
thin (1-4 particle diameters thick) hill-type structures (fig. 8-3, b). All 
motion at these frequencies occurs by particles skirmishing over the top of 
such hills; the area between these thin hills remains clean. This particle 
formation forms to minimize the reluctance in the region below the field 
structure. The dominant mode of operation at lower frequencies involves walking 




Instantaneous Field at Point P With Time 
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Observer at P sees counter clockwise rotation 


Figure 8-1 Circulating Field above Traveling Wave Structure 


Travelling Wave Direction 



Figure 8-2 Instantaneous Picture of Agglomerates Under the Travelling Wave 
Structure in the Duct 


802 - 


-209- 

Dlrectlon of Travelling wave 

< 




^ Time 


(a) Single Agglomerate Walking Motion (f<40 Hz) 



(b) Stationary Particle "hills" Formad above 60 Hz 


Figure 8-3 Particle Formations 




...i-.- ..-is*.- ■*»-. '.n.' • .., ■ 


- 210 - 

unifoTmly on the duct wall with a wave undulation moving opposite to the particle 
migration in the same direction as the travelling wave. 

A 

fin alternative step for studying these effects is shown in fig. C8>4) . 

A 4-pole, 3 phase synchronous stator now serves as the field structure. A 
cylindrical chamber with a rectangular plate of plexiglass suspended on pins 
serves as the secondary member. It is sealed and filled with water and a small 
amount of iron powder. A clockwise field causes the powder to circulate counter 
clockwise- on the walls of the chamber. The powder drives the water by viscosity 
vhich in turn causes a rotation of the paddle. The experiment is especially 
suited to correlating particle speed with frequency. The speed was found to 
increase somewhat linearly with frequency up to 20 Hz (indicating that 
agglomerate length is a nonlinear inverse function of frequency, decreasing 
sharply after 20 Hz) . 

These effects indicate particle ducting via a travelling wave structure 
is feasible. The exact nature of drag reduction benefits obtainable from 
particles confined and shuffled in this manner is not clear. The particles 
do bounce and rotate along the surface and provide a mechanism for both trans- 
fering momentum from the wall to the flow or vice versa. This transfer can 
be coherent in the sense of tending to prolong a net circulation, but it 
appears that it could also function on the scale of the turbulent eddies. A 
longer interaction in the free stream flow might be obtainable by actually 
releasing the particles by using a standing, rather than travelling wave field. 
The consideration of particle loss may cause this mode of operation to be im- 
practical. In any event, operation in either mode must be at a frequency less 


A 

This experimental apparatus was built and studied by Ed Wooten, Massachusetts 
Institute of Technology. 
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than 60 Hz. This restriction necessitates the use of a long pole pitch if a 
travelling wave is used to synchronously (a distance 1/3 of the way into the 
boundary layer where the speed is roughly half the free stream speed) shuttle 
particles along the aircraft skin. A 500 MPH jet operating at 40 Hz would 
require a wave structure pole pitch 


P 


• h 


(500) (-5280) 



1.4 m 


(3-1) 


The results of this thesis are directly applicable to the design of a precipitator . 
for capturing the particulate at the aircraft tail and shuttling it to the nose. 

As energy needs ejqiand, it may be desirable to apply this process to ships 
and underwater vessels. The travelling wave interaction could in fact be the 
basis for propulsion » perhaps with a reduction in hydrodynamic noise associated 
with conventional propulsion. The restriction on weight would not be as severe. 

The smaller difference in density between particulate and flow medium will be an 
advantage in terras of particle loss. The implantation of a travelling wave 
structure near the skin of the vessel should be an easier task as well. 

C. Synopsis and Related Research 

The thrust of this thesis centered on the prediction of magnetic particle 

t 

ducting and precipitation in turbulent air flows. Towards this end two models 
were developed, one appropriate for light particles (< 20u) and the other for 
heavy particles (> 20y) . The light particle model requires the solution of a 
diffusion equation with appropriate boundary conditions ijiposed where trajectories 
enter the duct volume. Inertial effects are unimportant. The heavy particle 
hybrid inertial-diffusion model represents the particle momentum balance in 
Lagrangian co-ordinates with an additional diffusion term added in by super- 
position. In both models, the effects of turbulent diffusion are lumped into a 
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measurable diffusivity constant. By means of correlation with experiments, 
both models have been shown to provide a useful degree of accuracy in predicting 
particle precipitation. Significant improvement over analytical models was 
shown in correlations with experimental data. Tliis encouraged the application 
of the hybrid diffusion model to the problem of particle containment in the 
aerodynamic boundary layer. Results indicated about 4% of the injected particulate 
would be lost if permanent magnet collection were used in a boundary layer of 
typical aircraft. The analysis further revealed 90% of the particulate would 
be precipitated over a 4 to 6 meter length of the aircract. Thus, unless it 
were reintrained in some manner of benefit to drag reduction, it also would be 
lost as a drag reducing agent. 

The analysis indicated that to improve correlation of experiment with 
predictions, the three dimensional nature of the problem must be brought into 
the model. Specifically either the flow and field structure needed to be made 
wider, or the three dimensional considerations needed to be added to the model. 
Incorporation of variations of flow over the width of the channel and viscoiis 
effects of the side walls into the model would lead to higher precipitation 
predictions in both heavy and light particle analyses. >lagnetic field edge 
effects vdiich are a function wavelength should explain some wavelength trends 
observed. For application to precipitation technology, for example coal 
desulfurization, the consideration of small gradients around the tips of 
precipitated dendrite structures may be the most significant area for further 
study. The dependence of precipitation on loading would give a clue as to the 
magnitude of this effect. 

The precipitation models should be of use to particle ducting studies in 
which loss of particulate is undesirable. The models developed should aid 
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researchers in particle augmented drag reduction since precipitation would 
noxmally be required sin5)ly to conserve particles. 

The wind tunnel is quite useful for experimentally studying general drag 
reduction effects and in this context, the above research is directly applicable 
in collecting particulate. 
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i\ppendix A * Particle Magnetizable Spherical Particle 

It is helpful in determining the nature of the magnetic force to consider 
the analogous polarized particle in an external electric field. The polar* 
ization originates from charges separated by a distance d (fig. A*l) , and the 
force can be examined separately on the two charges, i.e., 

1 - 2E(r + a) - 2E(r) (A-1) 

En^loying Taylor's expansion and assuming d small, (A*l) becomes 

r- 2(FCr) + ar.vE) - 2E*(r) (A-2) 

or 

7 » p-vE ^ (A-3) 

Where p is the particle polarization given by 2^. The material can be 
considered to have a number of dipoles per unit volume n^ . From equation (A-3) , 
the force density F is had by multiplying by np, i.e., 

F = n^ -vF = F-vF (A-4) 

The force on a magnetized particle now follows by analogy. No magnetic 
raonopoles exist, but one can certainly examine magnetic dipoles and define M 
as the number of magnetic dipoles m per unit volume . Now the electric force 
of (A-4) can be obtained via energy arguments, i.e., by taking the gradient 
of energy of electric dipoles in an external electric field. It is legitimate 
to use the same energy arguments with magnetic dipoles and the result must be 
the same, i.e. , 

F = (A- 5) 

where the u comes in because of the historical definition of M 
0 

(V#yQ(5T + H) = 0 *.v-(EqF + F) =a)- ^ 

The requirement of a gradient in the external field is intuitively reasonable; 
there must exist a different field at either end of a particle (acting 
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differently on the poles) to yield a net force. 

To apply (A-S) in detexmining the force, it is necessary to examine the 
field internal to a highly permeable sphere in an external field FT^ (fig. A- 2). 
For the external Z directed field FT^ and a particle of radius a, the scalar 


potential outside the sphere is represented as 

'Pq ■ -H^r cose + ^ cose 
r 

(A-6) 

'P.^ • Ar cose 
in 

(A- 7) 

viiere FT ■ -VH' 



Matching tangential H and normal B at r > a gives 


D - a* 


i(tj 

e-y 
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In a linear permeable material, the magnetization M and FT are related as 
Using (2-10) tlien, one finds M to be 


CA-9) 


M 


2 +H- 


a-.. 1 


FT. 


(A-10) 


With (2-5) and the fact that the particle force is obtained by multiplying by 
the volume, one finds that 
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or 


F 




j 7(R*H) 


F - a’7 (FT.H) 


Here the R due to magnetization has been brought inside the gradient, 
should be emphasized that the FT in (A- 11) is the external field only. 
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Appendix B - Refinement of Basic Precipitation Models 

The refined models and laminar flow test calculation are developed from 
the results of chapter. 2, section C. The -easiest model to refine is the Deutch 
model for a turbulent, fully mixed flow. The model posed thus far incorporated 
migration due to magnetic force only. Adding gravitational migration for 
particles with mass m ■ ' ters (2-21) to 



The solution of (C-1) again gives an exponential fall for the density. Since 
precipitation occurs at y - o, the density becomes 



x-x^ 

0 

2a U(mg gk) 
B 


) 


(B-2) 


Two additions to the particle trajectory laminar flow are in order. The 
first is the inclusion of the gravitational migration. The second deals with 
the replacement of the assumed uniform air velocity plug flow with the 


characteristic turbulent flow profile (U = i^) . The net particle fiux 

now contains a more accurate convective term along with two migration terms. 

eOce-^*^ 
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Following the procedure above of taking the divergence of f gives 
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(B-4) 


or in the particle's frame, 
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dn 

3F 


-a* 


g— * 


n 


along 


1 



As in the previous study it is noted that if the force field is not divergence 
free, (V*(VlT*H) ^ 0), the density is not constant along a trajectory line. 

For a single harmonic sinusoidal field (C-5) becomes 


dn , 

g 


n 


along 



ake-2ky 


(B-6) 


The particle trajectory dependence is (C-5, b) as before. The vertical tra- 
jectory's dependence on time can be found by integration immediately as 


mg 2kmg 


\mg + oke’^^o/ 


t 


(B-7) 


It is not as fruitful to seek the trajectory and density time dependence in this 
refined model, but is more convenient to solve for the x - y trajectory 
character and then to find the density dependence on a space dimension. Towards 
this end, dt can be eliminated from the two equations in (2-25, b) to give 

1 

U 1 ^ 

dx ® 

By mg ak _-2ky 
6 5 “ 


(B-8) 
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The density y dependence is obtained by substituting for dt its 
equivalent in terms of dy from (2,25»b), i.e., 


dt ^ 



and thus 


dn ^ 2ak^e~^^n dy 
^ iDg ♦ oke-2'^ 

It follows from (2-41) that 


(B-9) 


(B-10) 


a_= In / mg ♦ 

"o \mg*cke-2^ I 
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For any starting position yo, the density upon impact at y»0 follows directly. 
The particle's trajectory follows by numerically integrating (C-9). The pro- 
gram KDTRAG (Appendix C) uses a forward Eulerian integration, stepping y from 
yo to the duct wall at y=0. Fig. (B-1) shows typical trajectories of iron 
particles where the magnetic migration is upward. The calculation for these 
profiles assumed the particle size and magnetic field were such that the 
magnetic force balanced gravitational force midway up the duct. The increased 
downward curvature of the gravitationally dominated trajector/- lines near y»0 
reflect the convective (^)7 dependence. The magnetically dominated trajectories 
exhibit an even greater curvature (upward) due to the exponential field depend- 


ence. 










Appendix C - KDDTRAG 
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rMoP r.s Po/,;; 


This program calculates particle trajectories according to the re- 
fined Inertla-free model developed In appendix B. X and Y represent 
the particle^ 8 axial and vertical position respectively, while n repra 
sents particle density. The program Is In basic. 


10 DIH YO<3)»X0<3)fN<3>»Y<3»100)»X(3»100)»J9<3> 

12 PRINT ■THIS PROGRAM COMPUTES INERTIAL FLOW TRAJECTORIES 

13 PRINT *FOR SINGLE HARMONIC B— *U0< 1/9) rB» AND GRAVITY • 

14 PRINT 'T-GRAUfT-INfMAG »T-RESr K FOLLOW 

15 INPUT GlfMlfRlfK \ PRINT GlfMlrRlfK 
20 D2«.l \ D3».01 

30 FOR lal TO 2 

35 Y0<1>»«4252 \ Y0(2)*«315 \ X<I»0)=0 

40 FOR J*0 TO 400 

50 Y<IrJ+l)«Y0<I)-J)KD3 

55 IF Y(I»J+l)<1.00000E-03 GO TO 85 

60 U=G1/Rl)«c<2#< 1-Y< I »J+1 >>>''. 11 1111 

65 IF YdrJ+lX.S THEN U»Gl/Rl*<2)|cY< I » J+1 ) )''♦ 1 11 1 11 

70 F*U/<1+G1/M1*EXP<-25|CK*Y<I»J+1) ) )*D3 

75 IF J*0 THEN F=F*.5 

(30 X<I» J+l)«X<If J)+F \ NEXT J 

85 N< I )* < 1+G1/M1*EXP( -2!(CK*Y0< I ) > >/< 1+Gl/Ml > 

88 J9<I)*J 

90 XO<I)=X(I»J) \ NEXT I 
93 PRINT 'THE YOrXOf N FOLLOW 
95 FOR 1*1 TO 2 

100 PRINT YO<I>fXO<I)?N(I)» \ NEXT I \ PRINT 
105 INPUT HI \ PRINT 'THE X»Y TRAJECTORIES FOLLOW 
110 FOR 1*1 TO 2 
120 FOR J»1 TO J9<I) 

125 PRINT X<I»J) JY(I»J)» \ NEXT J \ PRINT 
130 INPUT H2 \ NEXT I 
135 STOP 
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)\ppendix D - Analytical Diffusion Analysis 

The point of this analysis is to determine the nature of the density 
gradient at a boundary separating regions at high and low diffusivity. Figure 
(D*l) shows the problem considered here, two thin layer regions of small 
diffusivity d and thickness 6 separated by a large core region of diffusivity 
D and thickness (2a~25). Only gravitational migration exists and no x depend- 
ence is considered. Injection of particulate occurs at the upper wall and the 
boundary conditions are n ■ n^ at y = 2a and n » 0 at y » 0 (perfectly absorbing 
lower wall) . 

The flux in region 1 is 


r - -d 
y ^ 




(D-1) 


In the steady state with * 0 » 
2 

3 n _ mg_ 3n _ 

V B 17"“ 


CD-2) 


The same equations apply to region 2 with d->D. Noimalizing to the duct width, 
we have 




(D-3) 
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Figure D-1 Analytical Diffusion Study Model 



CD- 2) becomes 
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The general solution of (D-4) is 


n * + C2 e 
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(D-S) 


With the boundary condition n = 0 at ^ » 0, we have 

-Id Z 

^g ’ 

~1 * ^1 ■ ® 


(D-6) 


Similarly for region z, the density is 

^D „ 


n2 » C3 + C4 e 




(D-7) 
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Imposing the condition of continuity of particle density and flux at = 


deteimines n, to be 
—z 


i\ 


n = c, 1 -le 
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Finally in region 3, we have a characteristic density 

I 


■g 


1I3 = C5 + cg e 


(D-9) 


and the boundary conditions n » 1 at 21 = 1 and continuity of density and flux 
at = (1-^ . Combining (D-6) , (D-8) , and (D-9) with these conditions yields 
the densities for each region as 


1 # 
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where 
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1 - e 


"g 


( 1 - 26 ; 


.-'d 


(D-IO) 




Figure (D-2) shows a typical density profile for this problem assuming 
D»d. The ratio of density gradient in the core region to the boundary layer 
region at ^ or 1-^ is found straight away from (D-10). The result is 


3n I 

^ region 2 

lE 

region 1, 
or 

region 3, }^l-5 


Id 

^d 


(D-ll) 


Thus, for small diffusivity d, the density gradient is nearly zero in the 

*^D 1 

core . This is explicityl shown for the two curves with ~ The deished 

Id - 1 ® 

curve with ” IM displays a sharp vertical slope at y=»S. Also the effect 
of smaller core diffusivity is seen to squeeze all the particulate into the 
lower layer. 
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Appendlx E <- KDDI4 ; KDD19 




heproducibility 

ORIGINAL PAGE IS 


OP THE 
POOR 


These fortran programs calculate particle density and deposition 
according to the approximate diffusion model presented In chapter 4-B. 
KDDI4 Imposes the boundary condition n*0 on the upper boundary , 
while KDDI9 uses a zero density gradient condition there. 


KDDI4 


CALCULATION OF IiUCT CONCENTRATION PROFILE 
SMALL PARTICLE TURBULENT DIFFUSION CASE 

DIMENSION B(10»80>»X(80>»AN(80) rBO< 16> »AM( 10) »U( 10) » AMK< 10) 

COMMON AK f BO » U1 f AMMO f NHARM » DLTUB 

URITE(Sr^) 

5 FORMAT THE NUMBER OF ROWS AND COLUMNS ARE') 

READ(Sf 10)IR0Uf ICOL 
10 FORMAT < 215) 

URITE<5»15)IR0UrIC0L 

15 F0RMAT(2I5»/» ' INPUT VARIABLES G1»AK»XB»AMMV»DLXMS»XMAG»DLTUB' ) 
READ ( 5 f 20 ) G1 » AK > XB » AMMV r DLXMS t XMAG i DLTUB 

20 F0RMAT(F15.2) 

URI TE < 5 > 25 ) G 1 f AK » XB » AMMV f DLXMS » XMAG f DLTUB 
25 FORMAT (' THE INPUTS G1 »AK»XB» AMMV r DLXMS » XMAG » DLTUB ARE' f /6F10.4) 
READ(5r23)AHl 
23 F0RMAT(F10.2) 

URITE(5f 16) 

16 FORMAT<' THE ♦ OF HARMONICS AND THE NORMALIZED B-FIELDS'/' 

1 FOLLOW— B < N ) s''H*J|t2*e5tcK/b/ < D/2a ) ' ) 

READ(5f 10)NHARM 
READ(5f20)(B0<I)fI=lf NHARM) 

WRITE <5r 17) <B0< I ) » 1=1 »NHARM) 

17 FORMAT <8F10. 5) 

DLX = 1./<IC0L+1.) 

W=12.7 

DLY = <1.-2.*DLTUB/W)/IR0W 
U1=.5968/XB 
YU«1.-DLY-DLTUB/W 
XPOS=XMAG + DLX*XB 
CALL DT2A(XP0S»DT) 

CALL U0AM<YU»U2»AM2»AMK2> 

U(1)=U2 

AM(1)=AM2 

AMK<1)=AMK2 
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G2-61/DT 

AM2-AH2/DT 

AMK2-AMK2/DT 

U2-U2/CIT 

B < 1 » 1 ) « 1 . -DLX/IJ2* < AMK2+ < 62+ AM2 ) /2/HL Y+ 1 . /DL Y#3lc2 > 

BO 30 I*2>1R0U 

YU»1 *-BLYJ|cl-DLTUB/W 

CALL U0AH(YU>U2»AM2f AMK2) 

U<I)«U2 

AM(I)-AN2 

AHK(Z)*AMK2 

AM2bAM2/DT 

AMK2-AMK2/BT 

U2-U2/BT 

B ( I » 1 ) » 1 . -0LX/U2HC AMK2 
30 CONTINUE 

BO 50 Jb 2»1C0L 
XPOSaXMAG + J«BLX«XB 
CALL BT2A<XF'0S»DT) 

BO 45 lair IRON 

G2a61/DT 

AM2»AM<I)/DT 

AHK2aAHK(I)/BT 

U2aU<I>/BT 

IF(I.GT.l*ANB.I.LT»IROW>B<l»J>aB(I»J-l) + BLX/U2)|c( 

1 -AMK2J|cB<I»J-l) + (62+AM2)*<B<I-l» J-1)-B<I + 1» J-l))/2/BLY 

2 +<B<I-l»J-l)+B<I + lf J-1>-2*B(I» J-1) )/BLY5 |c#2) 

IF<I*EQ.l) B<I»J>aB<I»J-l) +BLX/U2*<“AMK2»B<If J-1) - 

1 ( Q2+ AM2 ) #B < I + 1 r J- 1 > /2/BL Y+ < B < I + 1 » J- 1 ) -2*B ■ X » J- 1 ) ) /DL Y**2 ) 

IF(I.EQ.IR0W)B<I»J)aB<I»J-l>+BLX/U2*<-AMK24cB<If J-i,> 

1 +<2>lcB<I-l» J-l) - 2.:lcB(I»J-l) >/BLY!«cj|c2) 

45 CONTINUE 
50 CONTINUE 

BO 111 Jai,lC0L»10 
URITE(5f?0}(I»Ial»10)f J 

90 F0RMAT<' THE DENSITIES ARE AS FOLLOWS" »/»3X» II f9<6Xr 12) » 

1 /»" THE FIRST COLUMN IS"»I3) 

JJaJ+9 

IF ( J J ♦ GT ♦ I COL ) J Ja I COL 
DO 108 lalfIROU 

WRI TE < 5 » 1 00 X B < I » J J J ) » J J J* J » J J ) 

100 F0RMAT<10F8.4) 

108 CONTINUE 

READ(5rll0)AH3 

110 FORMAT (FIO. 2) 

111 CONTINUE 

YU a DLTUB/W 
DLYO a YU/10. 

CALL U0AM(YU»U2»AM2»AMK2) 

AM0aAM2 

SUMa«5*U2/(Gl+AM2) 

DO 115 J=lf9 
YUaYU-DLYO 

CALL U0AM(YUfU2f AM2»AMK2) 
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SUM-SUN-fU2/ ( G 1 + AM2 > 
IIS CONTINUE 

SUN-SUN»OLYO 

URITE(S>110)8UM 




CALL UOAN ( 0 • » U2 » AN2 > AMK2 ) 

DO 120 J>lflCOL 

AN(J)«G(IRONr J>)K(01+AN0)/(GlfAM2>«DLXHS 
X<J>-(DLX«J<fSUH)«XD 
120 CONTINUE 

URITE(S»125>(X(I)»AN(I)»I«lrIC0Lr3> 

125 F0RMAT(4E15*4> 

STOP 

END 

SUBROUTINE UOAN( YU»U»AM» ANK) 

DIMENSION BO(U) 

COMMON AKrBOrUlfAMMUrNHARMrDLTUB 


Y»YU 

IF<YU*GT« .5)Y»1.-Y 

U-U1«(2*9|cY>«:K*1111111 

AM«0« 

DO 15 I>l»NHARM 
P»2*3|CI3|CAK*YU 
IF<P.GT.20*>G0 TO 15 
AM»AM+AMMV*BO < I ) I #EXP ( -P ) 

15 CONTINUE 

IF<YU.LT.DLTUB/12.7>G0 TO 30 
AMKsO* 

DO 25 I«lrNHARM 
P«2*#rtcAK»YU 
IF<P»GT*20. )G0 TO 25 

AMKsAMK+2 * « I ««2«AMMV«B0 < I > ««2« AK«EXP ( -P > 

25 CONTINUE 
30 RETURN 
END 

SUBROUTINE DT2A<X»DT) 

DX»*1524 

IF<X«LT.DX)DTa *003625 
XEls2.«0X 

IF < X ♦ GE . DX . AND . X ♦ LT . XE 1 ) IiT« . 003625+ . 020915* < X-DX ) /DX 
IF<X.GE.1.33985)DT=. 00178 
IF(X«LT.XE1«0R*X«GE*1. 33985) 60 TO 50 
~ DT*. 054658“. 15085ntcX+.2<Xr784*X**2-*084262#X**3"* 07902* 
1 X**4+ . 028938*X**5+ . 065659*X**6- . 033616*X**7 


50 DT«DT/*127 

RETURN 
END 


n nnnci 


KDDI9 
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CALCULATION OF DUCT CONCENTRATION PROFILE 
SHALL PARTICLE TURBULENT DIFFUSION CASE 
CHANGE UPPER BOUNDARY CONDITION 

DIMENSION B<10»80>rX(80>rAN(80)fB0(16)»AH<10>rU(10) rAHK(lO) 

COMMON AKfBO»UlfAMMV»NHARM 
URITE(5f5> 

5 FORMAT (' THE NUMBER OF ROUS AND COLUMNS AREM 
READ<5»10)IR0U»IC0L 
10 FORMAT (215) 

URITE(5f 15)IR0Uf ICOL 

15 F0RMAT(215»/i ' THE INPUT VARIABLES 61 t AKf XBf AMMVf DLXMSf XMAG ARE') 
READ ( 5 f 20 > 6 1 » AK r XB f AMMV » DLXMS $ XMA6 

20 F0RMAT(F15*2) 

WRITE ( 5 f 25 ) 61 r AK f XB f AMMV » DLXMS » XMA6 
25 FORMAT (' THE INPUTS 61 »AK»X6f AMMV» DLXMS >XMA6 ARE' f /6F10*4> 
READ(5»23)AH1 
23 F0RMAT(F10.2) 

UR1TE(5»16> 

16 FORMAT(' THE ♦ OF HARMONICS AND THE NORMALIZED B-FIELDS'/' 

1 FOLLOW~B ( N ) ■"H##2*eJlcK/b/ < D/23 ) ' ) 

READ(5tlO)NHARM 
READ(5»20)(B0(I)f 1»1»NHARM> 

URITE(5»17)(B0(I)»I»1»NHARM> 

17 FORMAT <8F10. 5) 

DLX s 1*/(IC0L4-1. ) 

W=12.7 

DLY * <1.~1./W)/<IR0U-1* ) 

U1»*5968/XB 
YU»l.-*5/W 
XP0S»XMA6 -f DLX«XB 
CALL DT2A(XP0S»DT) 

CALL U0AM(YU»U2>AM2»AMK2) 

U(l)sU2 

AM(1)sAM2 

AMK(1)bAMK2 

62*6 1/DT 

AM2*AM2/DT 

AMK2*AMK2/DT 

U2*U2/DT 

B < 1 » 1 ) * 1 . -DLX/U2* ( AMK2 ) 

DO 30 I *2 f IRON 
YU*l«-DLYHc<I-l. >-.5/W 
CALL U0AM(YUfU2»AM2f AMK2) 

U(I)*U2 

AM(I)*AM2 

AMK(I)*AMK2 

AM2*AM2/DT 

AMK2*AMK2/DT 
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U2-U2/DT 

B ( I » 1 ) -1 . ~BLX/U2«AHK2 
30 


G2-61/DT 

AH2«AN(I)/IlT 

AMK2»AMK<I)/DT 

U2«U<I)/HT 

IF<1.GT.1*ANIM.LT*1R0U)B(I» J)«B(I»J-1) f £iLX/U2«( 

1 -AMK2»B(It J“l)+<G2+AM2^*<B<I-lf J-1 >)/2/DLY 

2 +<B<I-1»J-1)+B<I + 1,J-1)-2*B<I» J-l>)/riLY!lc)|c2) 

IF<1.EQ«1) B(If J>«B<I»J-1> +BLX/U2#<-AMK2*B<I»J-1) 

1 +C2#B(I+lf J-l)-2.JlcB(If J-l> J/DLY*»2) 

IF < I ♦ EG . IROU ) B ( I » J > »B < I » J- 1 ) +DLX/U2)tc < - AMK2HCB < I f J- 1 ) 

1 +<2*B<I-1»J-1) - 2»#B<I»J-1))/DLY**2> 

45 CONTINUE 
50 CONTINUE 

DO 111 J«1»IC0L»10 
URITE(5»90XI»I»lf 10)»J 

90 FORMAT THE DENSITIES ARE AS FOLLOWS' »/f3Xf II i9<6X» I2> » 
1 /»' THE FIRST COLUMN IS'rI3) 

JJ«J+9 

IF(JJ.GT,ICOL)JJ*ICOL 
DO 108 I«1»IR0W 

WRITE<5»100)<B(I»JJJ>»JJJ=Jf JJ) 

100 FORMAT <10F8. 4 > 

108 CONTINUE 

READ<5f 110)AH3 

110 F0RMAT(F10»2) 

111 CONTINUE 
YU « .5/W 
DLYO a YU/10. 

CALL U0AM(YUfU2»AM2>AMK2) 

AN0»AM2 

SUM«.5JtcU2/<Gl+AM2) 

DO 115 Jal»9 
YU*YU-DLY0 

CALL U0AM(YUfU2rAM2rAMK2> 

SUMaSUM+U2/ ( 6 X AM2 > 

115 CONTINUE 

SUM=SUM)lcDLYO 

WRITE(5rllO)SUM 

CALL UOAM < 0 . » U2 » AM2 r AMN2 ) 

DO 120 JalrICOL 

AN(J)aB(IROWf J)!»c(Gl+AMO)/<Gl+AM2)5»cDLXMS 
X<J)a<DLX#J+SUM)!|{XB 
120 CONTINUE 

URITE(5f 125)(X(I>fAN(I)»I=lf IC0Lr3) 

125 F0RMAT(4E15.4) 

STOP 

END 


CONTINUE 
DO 50 J«2»IC0L 
XP0S«XMA0 ••• J«DLX«XB 
CALL DT2A<XP0S»DT) 

DO 45 I«1»IR0U 


tt OF 



15 


25 

30 



SUBROUTINE UOAH( YUf U» AHr ANK) 
DIMENSION B0<16> 

COMMON AKrBOvUlrAMMUfNHARM 
Y-yu 

IF<YU.eT..5)Y«l.-Y 

U-Ui«(2*«Y)He«*nillll 

AM»0. 

DO 15 I-1>NHARM 

P-2«»I»AK»YU 

1F(P*6T*20*)60 TO 15 

AM- AM-I- AMMV»BO < I ) iic«2« 1 «EXP ( -P > 

CONTINUE 

IF(YU.LT. .5/12*7)60 TO 30 
AMK-0. 



DO 25 I-lfNHARM 
P-2*«I«AK»YU 
1F(P.6T.20. )G0 TO 25 

AMK-AMK+2 . »1 :K«2:lcAMNV)|eB0 ( 1 ) 4C)|c2«AK:KEXP ( -P ) 

CONTINUE 

RETURN 

END 

SUBROUTINE DT2A<X»DT) 

DX-.1524 

IF < X ♦ LT . DX > DT- ♦ 003625 
XE1-2.«DX 

IF < X ♦ GE . DX . AND ♦ X . LT . XEl ) DT- . 003625+ . 02091 5* < X-DX ) /EIX 

lF(X.GE«l*33985)DT-.0017a 

IF(X.LT.XE1. OR. X.OE. 1*33985) 60 TO SO 

DT- . 054650- ♦ 1 50851*X+ ♦ 201 784#X*JK2- . 084262*XjK)|c3- . 07902)K 

X««4+ * 028938«X«9K5+ . 065659«X««6- . 03361 6«X»«7 

DT-DT/.127 

RETURN 

END 
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App«ndix F- KDDF2 

This fortran program calculates particle densities and particulate 
deposition according to the refined small particle diffusion analysis 
dlsscused in chapter 4-B. The exact dlffusivlty X-Y dependence » vertical 
densities are solved simultaneously, and space derivatives consistent 
with causality are used. -,rrv OF 

Aww — ' . 1 » i I 

KDDF2 




C calculation 0“ DUCT COMCFNT^ATlON OPOPI'.E 
C S^^ALL PARTICLE TUOflULENT 0I®TUSI0M CASE 
c Boundary condition*? conststant with causality 
C ZERO OEMSITY 5PA0IENT IHDQSED ON LOWER SURFACE 

r 

DIMENSION S(6Si fX («0) ,AN(0O) ,R0 (9) ,R2<65) fS?K(6S) •Y(6‘?i . 
ll)R(65) .AMO (65) ,RR(65,ftO) tACBS*?) 

DOUBLE PRECISION A,R 
COMMON AK*RO*W.AMMVfNMARM 
WPITE(6,5) 

q format (• THE NIJMRFR of ROWS ANO COLUMNS ARr») 
oEAO(5»lO) TROWflCOL 
To rORMAT(2lS) 

WRITE(6,15) IR0W,TC0L 

TS FORMAT (?I5«/* » INPUT T(?« AMMV* WV, XR« OLXMS aREM 
READ (5*20) TGfAMMV*WV,XROLXMS 
PO F0RMAT(5Fin,2) 

write (6,25) TG.XB,Ammv,WV*DLXMS 
05 FORmAT(* the inputs T0,X3, AMMV,WV,0LX^S ARe • , /6FlO . A) 
NHARMsR.OI 

IF(WV.EQ.R..0R.WY.F0.12.)G0 to 30 

AK=15.71 

XMaO=,4B9 

RO (1) =.080356 

RO (?) =.0027 

n0(3)=. 00392 

RO (4) =.00151 

B0(5)=.001 

RO (6) =.0000931 

90 (7) =.0007905 

P0(R)=. 000537^ 

P0(9)=. 000352 
00 TO 40 

30 IF(WV.E0.1?.)G0 TO 35 
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AKsq,97S 
XMAG=.3175 
RO (1) =.099763 
B0(2)=, 00193 
R0(3)=. 00497 
eO(4)=. 0009366 
50 (S) =.0015 
RO(6)=.00067B 
B0(7)=.001?6 
00(8)=. 0005 
00(9)=. 00095 
GO TO 40 
35 AK=6.650 

XMAG=.076? 

00(1)=. 08993 
00(2)=. 00156 
00(3)=. 025 
B0(4)=. 000267 
00(5)5.00128 
00(6)5.0000533 
00(7)5.00364 
00(0)=. 000267 
00(9)5.0001064 
40 OLX 5 1,/(IC0U 
WW=12.7 
W=.127 

HLY 5 (l.-,031/WW)/(TO0W) 
DO 45 I=l,TPOW 
AN0(I)=1. 


REPRODUCIBILITY OF THE 
ORIGINAL PAGE IS KK'K 


Y(I) si.-DLY*! 

45 CONTiNUI- 

00 50 l5l,TPOW 

CALL U0AM(y(I) .'am.aMK.U) 

02(1) =AM 
02K(I)5flMK 
U0(I)=U/XB 
50 CONTINUE 

C INPUT FINISHED— SFT UP FINITE ELEMENT rOJaTlONS 
DO 151 NCOL51.ICOL 
XPOS=XMaG ♦ DLX«xfl*NCOL 
OY50LY 

DO 100 N0OW=1,TROW 

call DT2A (X®OS.DT.Y (MoOW) fOTMAX) 

TI 050 T/W**? 


C 50 , 

IF(Y(NP0W) .LT..125)C=1 . 

IF(Y(NR0W) .GT.,875)C=-1. 

GPDT 5 C*0T'^AX/W»*?*TD*12. 7/(1. 5875-. 031)/2./DY 
A (MP0W»2) =TG»TTD«2/0y»*2* ( 1 . ♦T5*02 ( NPOW) /W) /Dy ♦ 
1TG«U9 (NROW) /OLX^TG»0?K’ (N20W) 

A (MPOW. 1 ) 5-TG»TIo/DY»« 2- ( 1 . ♦TG*P2 ( NJROW) /W) /f)Y-GODT 
A (MOOW.3) 5-TG»TTD/DY*«2* GROT 
0 (MOOW) =TG»UR (MRDW) /ni X«AN0 (VIOOW) 
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C 


loo CON TXNUS ^ u /\n)/ 1 >Y 

AUtUsO. 

AdPOWfU s A(TOOW,1) -T3*Tn/nv»*? -G^^OT 
A(IPOW*3)sO. 

single column of flements now peaoy fop TNvFPSTOM 
13P call lNVRSF(A,q,IPOW) 

141 AN(NGOL)sAMf>(TPOW)«OLYMS 
00 145 JJsl.TPOW 


REPRODUCIBILITY OP THE 
ORICIVAL TS Poor 


ANO(JJ) a B(JJ) 

BB{JJ*NCOL)sANn< JJ) 

145 CONTINUE 

X (NC0L)=NC0L*ni X*XP 
?51 CONTINUE 

C DENSITIES are now KNOWN - OJTP'JT FOLLOWS 
IF(TG.ME.?4.4«4) go TO 108 
WPlTE(6*9n) (Idsll^Tl ,20) 

F0RMAT(» TME OENSITIEP APE AS FOLLOWS » ,/ , 1 OX , T 2, 3 ( f X , T ?) ) 
DO loa Isl,lROW 

WPITE(6,107) Y(I) ,(SR(T,J) ,J=11,71,20) 

F0RMAT(F5,3,4F8.4) 

CONTINUE 
YU=,031/WW 
DLY0=YU/S. 


90 


107 

10 « 


call U0AM( Y'J,AM,AMK,U) 

AMOsAM 

SUM=,5»U/(W/TO^AM) 
po 165 J®1,4 
YUsYU-OLYO 

call U0AM(yu,AM, aMk.U) 

SUMsSUM4.U/(W/TG^AM) 

T65 CONTINUE 

call U0AM(0, ,AM, AMKfiM 
PATs (W/TG+AMO) /(W/TG+ aM) 

SUM = SUM »0LY0*W 
WPITE(6,187) SIJM.PAT 

T87 F0RmAT(* The X integration and N patio ARE».?F12,5) 
00 170 Jal,ICOL 
AN(J)sAN(J)»RAT 


X(J)=X(J)+SJM 
170 CONTINUE 

WRITE (6, 125) (Xd) ,aN(T) d = l,IC0L,3) 
T25 F0PmAT(4E15,4) 

STOP 

END 

subroutine UOAM(yu,AH,AM<,iJ) 

dimension aO{o) 

common AK,S0,W, AMMV,NMARH 

Y=YU 

IF(YU.GT.,5) Y=U-Y 
U=4,572*(2.*Y)*«. 11 Ilf 1 1 
AM=0. 


15 


*^0 
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DO 15 laUMHAPM 
Ps?,*I*AK*yU 
IF(P.GT.20,)Gn TO l«5* 

AMsAM4AMMV»aO (T)<»*2*I*eXP(-0) 
AMk=AmK*2.»I**?*4MMV»O0 ( I) *»2*AK/W»EXP(-P) 
CONTINUE 

IFCTU.GT., 0309/1?. 7) GO TO 30 
U » .242*«?»YM»,127/.000015 
PETUPN 
END 

SUPROuTINE OT?a(X,OT,v#DTMAX) 


jn®. 


OO'P' 


DXs.1524 

IF(X.LT,OX)OTs. 003625 
XE1*2,*0X 

lF<X.GE.Dx.4Nn,X.LT.XFl) OT = . 0036 ?« 5 *, 02 o 915 «»(X-DX) /OX 
IF(X.GE.1.339A5)OT=,00175 
IF(X. LT.XEI. OR. X.GF. 1 . 33985 ) GO TO 50 
0T=. 054658-. 15085 1*X*. 201 78 A»X**?-,09^?62«X**3-. 0790?* 
1 X<»*4*. 028938* X** 5 +. 065659 *X** 6-.0 336 16*X** 7 
50 DTMAX = OT 

IF(Y.LT.. 1 ? 5 )DT= 0 T*Y/. 1?5 
IF(Y.LE. . 030905 / 12 . 7 ) DT = o. 
lF(Y.GT.. 875 )Ot=OT*n.-Y) /. 1?5 

return 


END 

subroutine rNVPSE(A.B.IROW) 

DIMENSION a ( 65 * 3 ) fP ( 65 ) * 3 ( 65 ) •U( 65 * 2 ) *Y( 65 ) 
DOUBLE PRECISION A*B*n,U*Y 
U(l*l) = A(l*?i 
U(l* 2 ) = A(l, 3 ) 

DO 30 I = ?*IROW 

U(I* 1 ) = A(T.?)-A(I- 1 , 3 )*A(I, 1 )/U(I- 1 * 1 ) 
U(I* 2 ) = A(I* 3 ) 

D(I) = A(I, 1 )/U(T- 1 * 1 ) 

30 CONTINUE 
Y(l) = B(l) 

DO 40 I s ?.IROW 

Yd) = R(I)- 0 (I)*Y(I- 1 ) 

40 CONTINUE 

8 (IROW)= Y(IROW) /U(IROW*l) 

DO 50 Js 2 *TR 0 W 


.I=IR 0 W* 1 -J 

Bd) = (Yd)-U(I,?)*R(I^ 1 ) )/U(T*l) 
50 CONTINUE 
RETURN 
FNO 

//GO.SYSIM DO * 

55 80 


P 4.684 14,063 8 . 1,5 1.1875 

/* 


/*EOJ ******** 


r> r» r> n o o 
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- REWOCf 

Appendix G - KDDFl OBIGTN^^^ 

This fortran program calculates light particle densities and particu- 
late deposition consistent with causality. Both the x and y diffuslvity 
dependence are included in the model. Diffuslvity is allowed to go to zero 
at the laminar boundary layer and thust no boundary condition is imposed 
at the lower core boundary. 

KDDFl 


r&LCULATlON Qc- OUCT CONCFNT^ATIOM PPOPT.E 

small papticl- tuopulent ot^’^usiom case 

BOUNOAPY CONOITIOMS CONSTSTANT with causality 
NO boundary condition on lower SUOFACE 

dimension (80) •AN(80) ,B0 (9) *02(^5) tS?K(65) ,Y(6‘^) • 

1UR(65) .AN0(6S) .BB(A5,R0) «A(65f'») 

DOUBLE PRECISION A,0 
CONMOM AK,Bn,N* ammv/.nmapn 
WRITE(6,5) 

q FORMAT (• THE NUM8FP Or ROWS ANO COLUMNS APr») 

PEAD(5*10) TPOW.'ICOL 
)0 F0PMAT(2I5) 

WPITE(6,1S) IPOW,TCOL 

7s FORMAT (2lS«/« » INPUT TO*AMM\/*WV*XB*DLXMS &RF») 
PEAD(5f?0) TG»AMMV»WV,yBOLXMS 
PO FORMAT (SFIO. 2) 

WRlTE(6f25 ) TG.XR, Ammv,WV*DLXMS 
PS FORMAT (• THE INPUTS TO,X3^ AMMV,WV»OLXMs ARri ,/^Fl0,4) 
NHARM=9,01 

TF(WV.E0.8..0R.WV.E0.12,)GO TO 80 
AK=15.71 
XMAG=.489 
B0(1)=. 080856 
80(2) =.0027 
e0(3)=. 00392 
R0(4)=, 00151 
R0(S)=.001 
RO (6)=. 0000931 
RO (7) =.0007905 
RO (8) =.000537^ 

BO (9)=. 000852 
GO TO 40 

lF(WV.E0.1P.)Gn TO 35 


80 


o o 


240 


I ^ 



I 


I 


AKsg.qTS 
xMAG=.3175 
R0(l>s*f)99763 
31 B0(?)a. 00193 
B0(3)=.004B7 
B0(4) s. 000536ft 
B0(5)=*0015 
B0(6) =.000678 
B0<7)=. 0017ft 
90(8)=. 0005 
90{9)=.00005 
GO TO 40 
95 fiKsf.650 

XMag=. 0762 
R0n)a. 08993 
B0(2V=. 00136 
90(3)=. 025 
80(4)3.000767 
90 (5)3.00179 
90(6)3.0000593 
90(7)3.00364 
BO (8) =.000767 
^ 90 (9) =.0001064 

40 DLX = l./(TCOL) 

WW=l2.7 

W=.127 

OLY = (l.-.031/WW)/(T70W) 

00 45 Isl.TBOW 
AN0(I)=1. 

Y(I)=1.-DLY»I 
45 COMTiNUE 

00 50 l3l,T90W 

CALL U0AM(Y(I).AM,AM»<,U) 

B2(I)=AM 
B2K(I)=AMK 
UR(I)=U/XB 
50 CONTINUE 

INPUT finished— SET UP FINITE ELEMFNT EQUATIONS 
KENT LOVES LANA 
00 151 ncol=utcol 
XPOS=XMaG ♦ 0LX»X8*NCnL 
DYaOLY 

DO 100 NR0i«f=l .TPOW 
CALL DT7A (XOOS.nT.Y(MOOW) »OTMAX) 

TI0=0T/W**7 

C* 0 . 

IF(Y(NR0W) ,GT..875)C=-1 . 

IF(Y(NR0W) .LT,.175)C=1. 

grot = C^»0T'^AX/W**?»TG*12. 7/(1,5975-. 031)/?, /ny 
A (NO0W.2) =TG»TT0»2/DY*-»2* ( 1 . ♦Tr,*®? (N^Ow) /W) /DY ♦ 
1TG«UR (NPOW) /0LX^TG»97 k (N^OW) 

A (NPOW. 1 ) =-TG*tio/OY«« 2-( 1 .*TG*92(NP0W) /Wl /OY-GPOT 
A (NPOW. 3) =-TG*TI0/nY*»2->3R0T 
9 (NPOw) =TG»UP (NPOW) /DI X*A^)0 (NPOW) 


Too CONTINUE 

A ( IROW# 1 ) sA { IROW* 1 ) ♦GOOT-OTMAX* ( .03l/WW*0V) /, 12‘5/OY*TG/W**?/nv*To, 
AdROW*?) sAdPOWf?) ♦OTMAX*(, 011 /WW* 0 Y) /, 1 ?5/0Y*T6/W**2/nY*l 0, 
Ad«l)sO. 

AdP0W.3)a0. 

C single column dp- elements NO<rf ready for INVEPSTON 
no CALL INvRSF(AtO.TROW) 

141 AN(NCOL)*AMO(TROW)*DLXMS 
00 145 JJ=1*IR0W 
ANOCJJ) s 9(JJ) 

9R(JJ*NC0L)=AM0(JJ) 

T45 continue 

X (NCOL) =NCDL*nLX»XR 
151 CONTINUE 

C DENSITIES ARE NOW KNOWN - O'JTPUT FOLLOWS 
lE(TG,ME.a4.6S4) GO TO 108 
WOITE(6*90) d,I=llf71,20) 

PO FOPMAT(» TME DENSITIES ARE AS FOLLOWS » ,/ • I OX , 12* 3 ( 6X , T?M 
DO 108 laldPOW 

WRITE(6*107) Y(T) * (RR(T*J) *Jall,7l,20) 

107 F0RMAT(E5*3f4FR.4) 

108 CONTINUE 
YUa,03l/WW 
DLY0=YU/5, 

CALL U0AM(YU,AM»AMK*II) 

AMOsAM 

SUMa,5<»U/('rf/TG^AM) 

DO 165 J=l,4 
YUaYU-OLYO 

call U0AM(yj*AM*AMK*U) 

SUM=SUM*U/(W/TG+AM) 

T65 CONTINUE 

CALL U0AM(0.fAM»AMK,in 
RATa(W/TG^AMO) /(W/TG+AM) 

SUM a SUM »0LY0«W 
WPITE(6*187)5Um,0AT 

187 FORMAT(» The X INTEGRATION AND N PATIO APE*.2E1P.5) 

DO 170 J=1,IC0L 
AN(J)=AN( J) *RAT 
X(J)=X(J)*SJM 

170 CONTINUE 

WRITE (6, 125) (X(I) *AN(T) d = l*IC0L*3) 

125 F0PmAT(4E15.4) 

STOP 

end 

subroutine UOAM(YU*AM,AM<fU) 
dimension R0(Q) 
common aK#RO»W.AMMV»NmaRM 
Y=YU 

IF(YU.GT..5) Yal.-Y 
(ja4.572*(2.*Y) **.11111 11 
AM = 0. 

DO 15 lal.MHAPM 
P=2.*I*AK*V'J 


-242- REJPRODUCIBILITy OP THE 

IF(P.GT.20.)GO TO 15 ORIGINAL PAGE IS POOR 

AMsAM*AMMV»eiO(t)**?»I*EX3(-0) 

AMKsAMK42,»I«h»?*AMMV*« 0 ( I)**?*AK/W*EyP{-P) 

Ts CONTINUE 

IE(YU.GT.. 0300/12, 7) GD TO 30 
U a .242 <h»?*YU*,127/. 000015 
30 PETURN 
ENO 

SUPPOUTINE OT2A<XtOT,YfOTMAX) 

0X5,1524 

IF (X.LT,DX)OTs, 003625 
XEls2,*DX 

lF(X.GE,OX.ANn.X.LT.XFl) DT s,003625**020915»(X-DX)/OX 

IF (x,GE.1 ,33955)075, 00175 

IF(X. LT.XEI. OP. X,GE. 1.33985) GO TO 50 

075,054658-, 150851*X*, 20 1794»X**?-,084262*X**3-, 07902* 
lX**4*.028938»X**5+.06e;659*X»*6-, 03361 6*X**7 
50 OTMAx 5 DT 

lF(Y.LT,,l?5)DTaOT*Y/.125 
IF(Y.LE,.031001/12.7) OT 5 o. 
lF(Y.GT..875)0Tsr)T*(l,-Y) /,125 

petuRn 

END 

SUBROUTINE INVPSF(A,B,IPDW) 

dimension a (65* 3) *B( 65) 0(65) *11(65*2) *Y (65) 

DOUBLE precision A*B*0*U*Y 
U(l*l) 5 A(l*2) 

U(l*2) 5 A(1*3) 

DO 30 I 5 P.Ipnw 

IJ(I*1) 5 A(I*2)-A(I-1,'3)»A(I,1)/u(I-1*1) 

IJ(I*2) 5 A(T*3j 

D(I) a A(I*l)/U(T-l*l) 

30 CONTINUE 
Yd) a'B(l) 

DO 40 I a ?*IPOW 
Yd) 5 Bd) -0(1) *Y(I-d 
40 CONTINUE 

P(IPOW)a Y(IROW) /U(IP0W*1) 

DO 50 J a ?,IP0W*1 
I 5 IP0W*1-J 

Bd) a (Y(T)-U(I*2)*B(I^1) )/U(T*l) 

50 CONTINUE 
RETURN 

end 

//go.sysin do 

55 80 

i?4.Q6 1.851P 12. 1,5 ,2389 

/* 

/*E0.) «•»*•»***•» 
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Appendlx H - KD1N4 


?%l8 fortran program calculates heavy particle Inertia trajectories 
and deposition according to the theory of chapter 4-C. No diffusion Is 
considered. Ntunerlcal Integration Is based on a forth order Runga-Cutta 
method. 

.TYPE KDIN4.F0R 
C 

C THIS PROGRAM CALCULATES THE TRAJECTORIES OF HEAVY 
C PARTICLES INERTIA DOMINATES— NO DIFFUSION 
C 

DIMENSION Y<3)-200> »YP(3r200> »X<3»200) ?XP<3»200) i-B0<15) » JMAX(3 

COMMON AMMQyNHARMf AK fBO 

WRITE(S»5) 

5 FORMAT<' GIVE ME T-ING T-INV ? T-RES ? AMMG ? YO » XPO » AK f DLT f DI A ' ) 
READ < 5 » 1 0 ) TG r T V » TR y AMMG y YO y XPO y AK y DLT y D I A 
10 FORMAT (FIS ♦ 5 ) 

WRITE ( 5 y 15 ) TG y TV y TR y AMMG y YO y XPO y AK y DLT y D I A 
15 FORMAT(' INPUTS TGyTVyTRy AMMGy YOyXPOy AKyDLTyDIA AREV8F10.4) 
WRITE(5y20) 

20 FORMAT (' GIVE ME THE tHARitlONICS AND THE B“S') 

READ(5y25)NHARM 
25 FORMAT (15) 

READ(5ylO) (BO(I) yI=lyNHARM) 

WRITE(5y30) <BO<I> yl^lyNHARM) 

30 FORMAT(' THE DENSITIES ARE ' y / y 8F10 . 4 ) 

GV*TG/TV 

GVR=TG*)K2/TR/TV 

Y(2yl)«Y0 

Y(lyl)*Y0+DIA/2. 

Y<3yl)=Y0“DIA/2. 

DO 35 I~ly3 
YP(Iyl)=0. 

XP<Iyl)*XPO 

X(lyl)=0. 

35 CONTINUE 

C INPUT IS FINISHEEiy VERTICAL INTEGRATION FOLLOWS 
DO 90 1=1 y3 
DO 85 J=2y200 

YO=Y(IyJ-l) - - - 

YPO=YP(Iy J~l) 

CALL AM(YOyAMG) 

DYP1=DLT>K(-GV>KYP0-1.-AMG) 

YHALF*Y0+DLT/2 . >KYP0 
CALL AM(YHALFyAMG) 

DYP2=DLT* < “GVjK ( YPO+DYPl/2 ♦ ) -1 . -AMG ) 
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YHALF= YO+DLT/2 . * < 2 . JKYP0+DYP2 ) /2 . 

CALL AM<YHALFi»AMG> 

D YP3aDLTJ« < -GV5K < YP0+DYP2/2 ♦ ) - 1 ♦ -AMG ) 

Y 1 *YO+DLT>k < 2 ♦ 3XYP0+D YP3 > /2 ♦ 

CALL AM(YlfAMG) 

n YP4aDLTJK < -GV* < YP0+DYP4 ) -1 ♦ -AMG ) 

YP < I f J ) -YPO+ < D YP 1 + 2 D YP2+D YP3 ) +D YP4 ) / 6 ♦ 

Y<I»J)«YO+ <YPO+YP<Ir J) )/2.>KDLT 
YHALF»<YO + Y<I»J))/2. 

C THE Y-INTEGRATIDN IS FINISHED r X-INTEGRATION FOLLOWS 
XPO-XPdr J-1) 

XO=X<Ir J-1) 

CALL UKYO»U»GVR) 

DXPl«DLT)K<-GV>KXPO + U) 

CALL UKYHALF»U»GVR) 

DXP2«DLT>K<-QV5i(<XPOfDXPl/2« ) + U) 

DXP3-DLTHC < -6V* < XP0+DXP2/2 ♦ ) +U ) 

CALL UKY(IyJ) fUrGVR) 

DXP4«DLT5K < -GV5K < XP0+DXP3 ) +U > 

XP < I r J ) aXPO+ ( DXPl +2 ♦ * ( DXP2+DXP3 ) +DXP4 ) /6 ♦ 
X<I»J)=XO+<XPO+XP<IyJ) )/2,3!<DLT 
IF<Y<Ir J) *LE»0* > GO TO 08 
85 CONTINUE 
88 JMAX<I)«J 
90 CONTINUE 

C X-INTEGRATION FINISHED r OUTPUT FOLLOWS 

DO no I=l»3 

WRITE<S?93>l!» JMAX<I) 

93 FORMAT Xi^Yi^UXyUY ARE FOR CASE- y JMAX— ' x 215 > 

WRITE<5y95) (Ji-Xa» J) i- Y ( I J ) ?XP <1 x J ) YP < I y J ) !» J-1 1- JMAX (I ) r 3.0 ) 
JM«JMAX<I> 

WRITE ( 5 X 95 ) JM X X < I X JM ) X Y ( I y JM ) y XP < I X JM ) x YP < I x JM ) 

95 F0RMAT<I5x4F15*5) 

READ(5xlO)AHl 
no CONTINUE 
STOP 
END 

SUBROUTINE AM<YxAMG) 

DIMENSION BO (15) 

COMMON AMMGxNHARMxAKxBO 
AMG«0. 

DO 10 I=lxNHARM 

P*2.)KAK5KY)KI 

IF<P*LT»0.)P=0» 

IF(P»GT»20) GO TO 20 
AMG=AM6+AMMG5KI>KBO < I ) >KN<2)KEXP ( -P ) 

10 CONTINUE 
20 RETURN 
END 

SUBROUTINE UKYxUxGyR) 

U=0» 

IF(Y»LT»0, )G0 TO 5 
IF < Y . LT . ♦ 5 ) U^GVR-'K < 2 * 'i< Y ) >!{>K . 1 1 II 1 1 1 
IF(Y,GEx ♦5)U~GVR>i<(2»>Ka ♦-¥) )>|oi<d.niin 
5 RETURN 
END 


Appendix I - KDIN7 
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This fortran program predicts particle densities and deposition 
according to the hybrld-lnertlal theory presented In chapter 4-D. 
Particle trajectories are calculated as In the diffusion-free model, 
but the diffusion term causes additional trajectory spreading. The 
analysis Is especially tailored to Intermediate sized particles where 
Inertia and diffusion are Important. 


CC 

C THIS PROGRAM CALCULATES THE TRAJECTORIES OF HEAVY 
C PARTICLES INERTIA DOMINATES— WITH DIFFUSION!! 

C 

DIMENSION Y(30)fYP<30)fX<30)»XP<30) rB0<9> » JMAX(30) » AN<30) » 

1 DLN(30> »AN00(30) »AN0(30> »6CLT(30> f GCLT2(30> r Y00(30) f X00(30) 

DOUBLE PRECISION Xi Yf XP»YP»AN»YOOf XOO»ANOO»XO»XPO» Y0» YPOr 
1 PCHY » PCHX f YCHAN » XCHAN » X3 r Y3 r RN3 
COMMON AMM6>NHARM»AK»B0 
URITE(5f5> 

5 FORMAT(' GIVE ME T-INGt T-INVf T-RESf AMMG» YOf XPOt WV»DLT' »/ 

1 ' DIA»XTUBE »BND-LAYER<CM)f tROWSf DIF-ON" > 

READ(5il0)TGfTV»TR»AMMGfY0>XP0»UVfDLTfDIA>XTUBErBL»R0UfD0N 
10 FORMAT <F15. 5 > 

URITE(5»15)TGfTVfTRf AMMGfYOrXPOrUVrDLTf DIAfXTUBE»BLrROUfDON 
IS FORMAT < ' TG » T V » TR » AMMG » YO f XPO i WV r DLT » DI A r XTUBE , BL » ROW i DON " 

1 f/>dF12.6) 

NR0U»R0U*f*01 
NHARM*9.+.01 
T*. 006244 

IF(WV.EQ*8* .0R.UV.EQ.12. >G0 TO 20 

AK»15.71 

BO(i)». 080356 

BC<2)*.0027 

B0<3)«* 00392 

B0(4>«. 00151 

B0<5)*.001 

B0<6>*. 0000931 

B0<7)=. 0007905 

B0< 8 >*.0005376 

B0< 9) *.000352 

GO TO 30 

20 IF(UV.E0*12. >G0 TO 25 

AK*9.975 
B0<1)*.099763 
B0<2>*. 00183 
B0<3>*. 00457 
B0< 4 >*.0005366 


B0(5}«*0015 
B0( 6) -*000678 
B0(7)-. 00126 
B0(8)-*0005 
B0(9>«*00095 
60 TO 30 

25 AK«6*649 
B0(1 >-.08993 
NHARM-l.f *01 
60 TO 30 

26 B0(2)-*001S6 
B0(3)-*02S 
B0(4>- *000267 
B0(5)-*00128 
BO (6) -*0000583 
B0(7)-*00364 
B0( 8) -*000267 
B0(9)». 0001064 

30 6V-T6/TV 

0VR«6V#T6/TR 
D6-T0*J(C2/* 127#3|c2/TV#riON 
0LY-0IA/(NR0U+1) 

BO 35 I-1»NR0U 

Y<I)«Y0+rHA/2*-lJ(criLY 

X(I)«0* 

XP<I)«XPO 

YP<I)=0. 

JMAX<I)-1000 

AN0(I)«<l.-ABS(Y<I)-Y0)/(DIA/2. > >XcjK. 1111111 
AN(I)-ANO(I> 

35 CONTINUE 

CALL DT(XTUBE»Ii) 

DLN < 1 ) — AN < 2 ) /2 * /DL Y/AN < 1 ) *D 

DLN < NROW ) -AN < NROW- 1 ) /2 . /DL Y/AN ( NROW ) #D 

DO 36 I-3»NR0U 

DLN(I-l) = <AN<I-2)“AN<I))/AN<I-l)/2./DLY!<cD 

36 CONTINUE 
DLXO-XPO*DLT 
INAX-NROU 

INPUT IS FINISHED* VERTICAL INTEGRATION FOLLOWS 
DO 120 J-lflOOO 
DO 40 I-lfIMAX 
YO«Y<I> 

YPO-YP(I) 

CALL AN(YOfANG) 

DYPl -DLT» ( -6VJKYP0- 1 . - AM6-D6#DLN < I ) > 

YHALF- YO+DLT/2 ♦ #YP0 
CALL AM ( YHALF »AM6) 

DYP2-DLT* ( -GV* ( YPO+D YPl /2 . ) - 1 . - AMG-DG)lcDLN < I > ) 
YHALF- YO+DLT/2 . # < 2 . *YP0+DYP2 ) /2 . 

CALL AM ( YHALF fAMG) 

DYP3=DLT;»c < -G VJ»c ( YP0+DYP2/2 ♦ ) -1 ♦ -AM6“DG:«cDLN < I ) ) 
Yl-YO+DLT# < 2 . * YP0+DYP3 ) /2 ♦ 

CALL AM<Y1»AMG) 
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DYP4«DLT* < -6V* < YP04DYP3 ) -1 ♦ -AM6-ri6#DLN < I > > 

YP ( I > ■ YPO+ ( DYPl +2 . ♦ < DYP2+DYP3 > 4DYP4 ) /A . 

Y<I)»YO+ <YP0+YP(I))/2.#DLT 
YHALF«<YO + Y<I>)/2* 

C THE Y-INTEORATION IS FINISHEDf X-INTEGRATION FOLLOWS 
XP0*XP<1) 

XO«X(I> 

CALL UKYOfUrOVR) 

DXPl«DLT*<-GV>lcXPO + U> 

CALL UKYHALFtUfGWR) 

»XP2«DLT#<-GV!|c(XP04riXPl/2« > + U> 

DXP3-DLT# ( -GV# ( XP0+DXP2/2 ♦ ) +U ) 

CALL Ul(Y(I)rU»6WR> 

»XP4»0LT» < -GV!K < XP0+DXP3 > +U ) 

XP < I ) «XPO+ < DXP 1 +2 . * < DXP2+HXP3 ) +DXP4 ) /6 . 

X < I > »XO+ < XPO+XP < I > ) /2 . *DLT 
IF<Y<I>.LE.O* > IMAX*IMAX-1 
IF<Y<I).LE.O. )JMAX(I)«J 
YOO<I>«YO 
XOO<I)»XO 
40 CONTINUE 

C X-INTEGRATION FINISHED »NOW GET DENSITY 
DO 90 I»1»NR0W 
XO«XOO<I) 

YOsYOOd) 

ANOO(l>sAN<I> 

IFd.EO.DGO TO 65 

IF<Y<I-1>»LE«0. .AND«JMAX<I).LT.J)GO TO 90 
IF(I.EQ«NROW)GO TO 68 
YCHANa<Y(I-l)-Y(I+l) >/2./DLY 
IF(Y(I+1>,LE*0*>YCHAN*(Y(I~1)-Y(I))/DLY 
XCHAN« < X < I ) -XO > /DLXO 
PCHY*(X(I-l>-X<I+l>>/2./DLY 
PCHX=»<Y<I)-YO)/DLXO 
GO TO 70 

65 YCHAN*<Y<1>-Y<2)>/DLY 

XCHANe < X < 1 ) -XO > /DLXO 
PCHY«<X(1)-X<2))/DLY 
PCHX*<Y(1>-Y0)/DLX0 

GO TO 70 - 

68 YCHAN»<Y<NROW-l)-y<NROW) )/DLY 
XCHAN= < X < NROW ) -XO ) /DLXO 
PCHYa<X(NROW-l)-X<NROW) )/DLY 
PCHXa < Y< NROW ) - YO ) /DLXO 
70 AN ( I ) s ANO < I ) / ( XCHAN# YCHAN-PCHY^ICPCHX ) 

90 CONTINUE 

C DENSITIES ARE KNOWN-GET D*6RADIENT<N) 

95 DO 110 1*1 » NROW 
DLN<I)*0. 

IF<I.6T.IMAX)60 TO 110 

IF<Y(I>.LE4BL/12«7«0R.Y(I>.GE. (1,-BL/12.7))G0 TO 110 
XPOS*X( I ))»:, 127+XTUBE 
CALL DT<XPOSfD) 

IF(I.EO.l) GO TO 100 
IF<I.EQ*NROW>GO TO 105 


RN3-AN(I^1> 

106 DLN<I)«<(AN<I-l)-RN3><c<X<I>-X00<I))-<AN<I)-AN00<I>># 

1 <X<I-1>-X3))/<<Y(I-1>-Y3)#<X<I)“X00<1)>-(Y<I>-Y00<I))* 

2 (X(I-1>-X3>) / AN(I)«D 
60 TO 109 

100 DLN<1)—AN<2)!K<Y<1)-Y<2>>/<<X<1)-X<2)>**2+(Y<1)-Y<2)> 

1 »»2)/AN(I>«D/2* 

GO TO 109 

1 05 DLN < NROW > -AN < NROW- 1 ) # < Y ( NROW- 1 > - Y ^ NROW ) > / ( ( X < NROW- 1 > ~ 

1 X<NROW)>#3|{2+<Y<NROW-l>-Y<NROW>))|c*2)/AN<NROW>#D/2« 

109 SIGN-1* 

ZF(IILN(I)*LT*0.)8IGN— 1* 

W-*127#S0RT<D/T) 

Z F ( S I GN«DLN ( I > * GT * U > ULN ( I > «S I GN«U 

110 CONTINUE 
ZF(Y(1>«LE.0.)60 TO 125 

114 IF(H0D(Jr5)«NE«0)00 TO 120 

112 URITE(5rlll)J 

111 FORMAT <' THE OUTPUTS-Xf Y»XPiYP» DENSITY iULN-CASE-' * 15) 

DO 115 I-1»NR0U 

URITE(5f 113>If)((I)fY(l)fXP(I)fYP(I) >AN(I) ^DLN(Z) 

113 F0RMAT(I5r6E12*5> 

115 CONTINUE 
120 CONTINUE 

C 

C NET INTEGRATION FINISHED--OUTPUT FOLLOWS 
C 

125 WRITE<5»126> 

126 FORMAT(' THE OUTPUT FOLLOWS V» 2X» ' I ' »X* 'JMAX' »6X» 'X' » IIX* 

1 'Y' f lOXf 'XP' »12X» 'YP' f 9Xr'GCLT' »SX» 'GCLT2' ) 

DO 135 I-lfNROW 

lF(I*NE*l*AND*l*NE«NROW>GDX-(X(I-l>>X(I-fl) )/2* 
1F(I.E0*1)6DX-X(1>-X(2) 

IF (I . EQ * NROW ) 6DX-X ( NROW- 1 ) -X ( NROW ) 

6CLT < I > -100 . / . 9/DI A#DL Y/GDXJtcWV/2 . /1 2 * 75«cANO < I ) 

6CLT2< I ) — 100 . *WV/2 . /12 . 7/ . 9/DI A/XPO^AN < I ) 3KYP < I > 

135 CONTINUE 

DO 140 1-1»NR0W 
X<I>-X<I)#.127 

WRITE<5»130)I»JMAX<I> jX<I)fY<I)»XP<I)»YP<I) r6CLT<I>*GCLT2<I) 
130 F0RMAT(2I4>6E12*5) 

140 CONTINUE 
STOP 
END 

SUBROUTINE AM(Y*AMG) 

DIMENSION B0(9) 

COMMON AMMG»NHARM*AKfBO 
AMG-0* 

DO 10 I-1»NHARM 
P*2.!|cAK5|cY*l 
IF(P.LT.O. )P*0. 
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10 

20 


s 


1 


50 


IF(P';0T.20.) GO TO 20 

AN6«AM6-f AMH6«1 !KBO ( I > )|c«2»EXP ( -P > 

CONTINUE 

RETURN 

END 

SUBROUTINE UKYrUfOVR) 

U«0* 

P-1./9. 




IF(Y*LE*0*>G0 TO 5 

IF<Y.GT*(l.-.5/12.7>>UsGVR»<2.)|c.5/12*7>!(t5|cp 
IF(U«NE*0.)G0 TO 5 


IF(Y.LT. *5>U-GUR»(2*«Y>«}|CP 

IF<Y.GE..5>U«GVR)|c<2.J(c<l.-Y))*)|cp 

CONTINUE 


RETURN 

END 

SUBROUTINE DT(XfD) 

Da « 061664- * 1 464 1 8«X-i- « 1 358S4»X»«2- * 01 3655)KXi|e«3 
- ♦ 03752#X*JK4- . 002257)lcXs|c:K5+ . 0099035ttX*)«c6 
RETURN 
END 



ui u ^ u ro f-j 

cn o ui o Cl (4 o cn M 
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Appendlx | - FOURIE 






This basic program performs a discrete fourier analysis of Che mag- 
netic field density waveform 1/4" above the permanent magnet structure. 
The analysis follows that described in chapter S. 


10 DIM X(3i)»B(16) 

PRINT "N HARMONICS FOLLOW— NEED 2N-1 DATA POINTS* 

PRINT *N«* \ INPUT N \ PRINT *X<I) FOLLOWS* 

P«2#N-1 

FOR I»1 TO P \ INPUT X<I) \ NEXT I 

PRINT *THE X'S ARE* \ FOR 1*1 TO P \ PRINT X<I)» \ NEXT I \ PRINT 

FOR J*1 TO N \ B(J>*0 

FOR 1*1 TO P 

B(J)*B(J>-i-X(I>/N)lcSlN(3*14159/N«J«I) 

NEXT I \ NEXT J 

FOR 1*1 TO N \ PRINT I»B<I) \ NEXT I 
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Appendix K " KD747G 


This fortran program calculates particle density and deposition 
for Intermediate sized particles injected over the skin of an aircraft 
flying 500 MPH at an altitude of 40,000 ft. The exact normal and axial 
dlffuslvlty dependence In the boundary layer Is used. Only magnetic mi- 
gration exists normal to the fluid stream flow; gravitational effects 
are Ignored, the analysis follows that described In chapter 7. 


CC 

C THIS PROGRAM CALCULATES THE TRAJECTORIES OF HEAVY 
C PARTICLES- INERTIA DOMINATES— WITH DIFFUSION!! 

C CALCULATES FLIGHT OF PARTICLES OVER 747 JET 
C ALTITUDE* 40»000 FT.» SPEED*500 MPH. 

C 

DIMENSION Y(30)»YP(30) »X(30) »XP<30)fB0<9>»JMAX(30>f AN(30) » 

1 DLN(30> » AN00<30) f AN0(30) rGCLT(30) »GCLT2(30) i Y00(30) f X00<30) 

DOUBLE PRECISION X»Yf XPt YPrANf YOOi.XOO» ANOO»XO»XPO> YO» YPO. 

1 PCHY f PCHX » YCHAN t XCHAN » X3 f Y3 f RN3 
COMMON AMMGrNHARM>AK»BO 
URITE(5fS) 

5 F0RMAT<' GIVE ME T-ING»T-INVf T-RES» AMMGi YO»XPO»WVf DLT' »/ 

1 ' DIArXTURBr *ROUS p DIF-ON ^ ) 

READ < S f 1 0 > TG f TV > TR r AMMG » YO » XPO f WV r DLT r DI A f XTURB » ROW f DON 
10 F0RMAT(F15.5) 

WRITE (5il5)TG»TV»TR»AMMGfY0fXP0fWVf DLT »DIA»XTURB»ROWf DON 
15 FORMAT< ' TGfTVf TR»AMMGf Y0»Af'0»WVf DLT»DIA»XTURB»ROW»DON^ 

1 »/»6F12.6) 

NR0W*R0W+«01 

NHARM*9.+.01 

IF(WV.EQ.8. .0R.WV.EQ.12. )60 TO 20 
AK*123.68 
B0<1 >*.080356 
B0<2)=.0027 
B0(3>*. 00392 
B0<4)*.00151 
B0<5)*.001 
B0<6>*. 0000931 
B0<7>*. 0007905 
B0< 8 >=.0005376 
BO <9 >*.000352 
GO TO 30 

IF<WV.EQ.12. >60 TO 25 
AK=78.54 
BOd > = .099763 
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B0(2)»*00183 
60(3>-*00457 
B0(4)-«000S366 
B0(S>«*0015 
B0<6)«» 000678 
B0(7)=*00126 
B0(8>».0005 
B0<9)**00095 
60 TO 30 

25 AK»52*36 

B0(l>».08993 
B0(2)s. 00156 
B0(3>»*025 
B0(4>« *000267 
B0(5)««00128 
B0(6>» *0000583 
B0(7)=*00364 
BO ( 8) s* 000267 
B0(9>s*0001064 

30 GVsi. 

6VR»TV/TR 
DG«TV#D0N 
DLYsDIA/<NROW+l) 

00 35 1-lrNROU 
Y< I >»Y0+DIA/2* -I*DLY 
X(I)*0* 

XP<I)=XPO 
YP(I)«0* 

JMAX(I)a850 
AN0<I)«<l*-ABS<Y<I)-Y0)/<DIA/2*> >3|t**lllllll 
AN(I)»ANO(I} 

YOO<I)=Y<I> 

XOO(I>*-XPO^DLT 

AN00(I)-AN(1) 

35 CONTINUE 

DLX0«XP0*DLT 

INAXaNROU 

INPUT IS FINISHEOf VERTICAL. INTEGRATION FOLLOWS 
00 120 JslfSSO 
00 40 I»lfIMAX 
YO=Y(I) 

YPO*YP<I) 

XO=X(I) 

XOT«XTURB+XO 

0L« . 37*XOT3|c < XOT/TR/ . 4701 )**(-. 2 ) 

GO TO 95 

39 CALL AM(YO»AMG> 

0YPl=DLT3|t < -GV)|cYPO“AMG-EiG)|criLN < I ) ) 
YHALF«Y0+DLT/2 . *YPO 
CALL AN(YHALF*ANG) 

0YP2=DLT* < -GV* ( YPO+DYPl /2 . ) -AMG-DG*OLN < I ) ) 
YHALF= YO+DLT/2 . * ( 2 . )|cYP0+DYP2 ) /2 ♦ 

CALL AM<YHALF»AMG) 




OP 
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DYP3*DLT* < - 6 V* < YP0+DYP2/2 . ) - AMG-D 6 *DLN < I > > 

Yl« YO+IiLT» ( 2 . #YP0+DYP3 ) /2 ♦ 

CALL AM<YlfAM 6 ) 

DYP4*DLT* ( - 6 V# ( YP0+DYP3 ) -AMG-D 6 «DLN ( I ) ) 
YP(I)»YP0+<DYPl+2*#<DYP2+DYP3)+DYP4)/6. 
Y<I)»YO+ <YPO+YP(I))/2,4cDLT 
YHALF«<YO + Y<I))/2« 

C THE Y-INTEGRATION IS FlNISHEIi» X-INTEGRATION FOLLOWS 
XPO«XP<I) 

CALL Ul<YO>U»GVR»DL) 

DXPl«DLT*<-GV)|cXPO + U) 

CALL UKYHALFfUfGMRrDL) 
DXP2=DLT»<-GV5|c<XPO+DXP1/2. ) + U) 

DXP3=DLT>|c < -GV# ( XP0+DXP2/2 4 ) +U ) 

CALL Ul(Y(I)rUfGVRfDL> 

DXP4aDLT5|c < -GV* < XPO+DXP3 ) +U > 

XP < I ) »XPO+ ( DXPl +2 4 * ( DXP2+DXP3 ) +DXP4 > / 6 ♦ 
X<I)bXO+<XPO+XP(I) )/24#DLT 
IF(Y<I)4LE40.> IMAXsIMAX-1 
IF<Y<I)4LE404>JMAX<I>«J 
YOO(I>*YO 
XOO(I)=XO 
40 CONTINUE 

C X-INTEGRATION FINISHED tNOU GET DENSITY 
DO 90 I«1»NR0U 
XO=XOO<I) 

YO«YOO<I> 

ANOO(I)-AN(I) 

IF<l4EG4l)G0 TO 65 

IF(Y<I-1)4LE4044AND4JMAX<I>4LT4J)G0 TO 90 
IF<l4EQ4NR0U)G0 TO 68 
YCHAN=<Y<I-1>-Y(I+1))/24/DLY 
IF<Y<I+1) 4 LE 4 O 4 )YCHAN«<Y(I-1)-Y<I) )/DLY 
XCHAN« < X < I ) -XO ) /DLXO 
PCHY=<X<I-1)-X(I+1> )/24/DLY 
PCHX=(Y<I)-YO)/DLXO 
60 TO 70 

65 YCHAN-<Y<1)-Y<2) >/DLY 

XCHAN* < X U ) -XO ) /DLXO 
PCHY=<X(1)-X<2> )/DLY 
PCHX* < Y (1 ) - YO ) /DLXO 

60 TO 70 ' 

68 YCHAN=(Y(NR0U-1)-Y<NR0W) )/DLY 
XCHAN= < X < NROW ) -XO ) /DLXO 
PCHY«(X(NR0W-1)-X<NR0W) )/DLY 
PCHX* ( Y < NROW ) - YO ) /DLXO 
70 AN< I) *ANO < I ) / ( XCHAN*YCHAN-PCHY!»:PCHX ) 

90 CONTINUE 

IF<Y<1) 4LE»04 >60 TO 125 
IF<M0D(Jf5) .NE40)G0 TO 120 
WRITE<5»91)J»IMAX 



non 


\r 
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91 FORMAK' Xfy»XP»YPfDENSITY»riLN“CASE«'fI 5 r ' 1 MAX«'» 15 > 

DO 93 l-lfNROU 

URITE( 5 f 92 >Z»X(I)»Y(I>»XP(l)rYP(I)rAN(I)>DLN(I) 

92 FORMATdSf 6012 * 5 ) 

93 CONTINUE 
120 CONTINUE 

C DENSITIES ARE KN 0 UN- 6 ET D» 6 RADIENT(N>— THE FOLLOWING STATEMENTS 
C ACT AS A SUBROUTINE TO GET THE GRADIENT 
95 DLN(I>aO* 

IF<I. 6 T.IMAX)G 0 TO 39 

IF<Y<I)* 6 E.DL« 0 R.Y<I)*LE* 0.)60 TO 39 ^,trmn.TTY OF THK 
CALL DT < DL » Y < I ) » D ) p AfiE IS ' 

IF<I.EQ* 1 > 60 TO 100 ORIGIN 

IF(I.EQ*NR 0 U )60 TO 105 
Y 3 »Y<I+ 1 ) 

X 3 aX(I+l) 

RN 3 aAN<I-M) 

106 DLN(I>a<<AN<I-l)-RN 3 >!|c(X<I)-XOO(I>)-<AN<I>-ANOO(I))* 

1 (X 00 <I- 1 )-X 3 ))/<<Y 00 <I- 1 )-Y 3 )#<X<I>-X 00 <I))-<Y<I)-Y 00 <I)) 

2 !|c<X 00 <I-l>-X 3 )) / AN(I)#D 
60 TO 109 

100 DLN(l)a-AN< 2 >)lc(Y<l)"-Y( 2 ) )/< <X< 1 )-X( 2 ) )*J|t 2 +< Y< 1 )-Y< 2 ) ) 

1 «« 2 )/AN(I)«D/ 2 * 

60 TO 109 

1 05 DLN < NROW ) =AN ( NROW- 1 ) * < YOO < NROW- 1 ) - Y ( NROW ) ) / ( ( XOO < NROW- 1 ) 

1 -X ( NROW ) > ** 2 + < YOO ( NROW- 1 > - Y ( NROW ) ) ** 2 ) /AN < NROW ) *D/ 2 . 

109 SIGNal. 

TaDL/* 7 / 223*5 

IF(DLN<I)*LT« 0 *)SI 6 Na-l* , 

WaSQRT<D/T) 

IF ( SI GN)lcDLN (I > . 6T * W ) DLN ( 1 ) =SI GN«W 
60 TO 39 


NET INTEGRATION FINISHED-OUTPUT FOLLOWS 


125 WRITE(5fl26) 

126 FORMATC' THE OUTPUT FOLLOWS V»2X» d ' »X» 'JMAX' »6X» 'X' » llXt 
1 "Y'flOXf 'XP'»12X»"YP"»9X»"6CLT'»8X» 'GCLT2') 

DO 135 lalfNRUW 

IF(I.NE*l.ANIM4NE*NR0W}6DXa(X(I-l)-X(I-l-l))/2. 

IF(I.EQ«l)6DXaX(l>-X<2> 

I F < I ♦ EQ . NROW ) 6DX=X ( NROW-1 ) -X < NROW ) 

6CLT < I ) *100 . / . 9/DI A*DLY/6DX*WV/2 * /12 ♦ 7*AN0 < I ) 

6CLT2 < I ) *- 1 00 . JlcWV/2 . / 1 2 ♦ 7/ . 9/DI A/XPO*AN < I ) # YP ( I ) 

135 CONTINUE 

DO 140 1*1 f NROW 

WRITE<5»130>I»JMAX(I)fX<I)fY(I)»XP<I) *YP<I)»6CLT<I)»GCLT2<I) 
130 F0RMAT(2I4f6E12.5) 

140 CONTINUE 
STOP 
END 
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SUBROUTINE AM(YfAN6> 

DZNENSION BO<9> 

COHflON AMMO»NHARH>AK»BO 
AMG«0* 

DO 10 Z»1»NHARM 
P«2.»AKS|CY«I 
IF(P.LT.O. )P«0. 

IF<P.6T*20.) 60 TO 20 
AM6>AM6-I-AMN6«I »B0 ( I > ««2«EXP ( -P > 

10 CONTINUE 
20 RETURN 
END 

SUBROUTINE U1 ( Y>U»6VR>DL> 

U«0* 

IF<Y.LE*0* >60 TO 5 
IF<Y»LT*DL)60 TO 3 
UsGUR 
GO TO 5 
3 P»l*/7* 

U«6VR*<Y/DL)**P 
5 CONTINUE 
RETURN 
END 

SUBROUTINE DT(DL»YiD) 

D=0. 

IF(Y.LE»0» .0R.Y.GT«DL)60 TO 50 
YN«Y/DL 

Da- ♦ 01031+2 ♦ 91226*YN+ ♦ 070557:f!YNJfc)|t2-13 ♦ 614712tYN*3ft3 
+4 . 179703»YN#*4+21 ♦ 233276*YNJlcJlc5-2 . 772269*Ym*6- 
15 . 559235* YN**7-8 . 465072* YN**8+ 12 . 027271 *YN**9 
DsD*. 0037*223. 5*DL 
IF(D.LT.O. )D=-D 
50 RETURN 
END 
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Appendix L - Refinement of Causal Fundamental Light Particle Model 

The density profiles predicted by the fundamental causal model are In- 
tegrated to give mass deposition per half wavelength on the lower duct sur- 
face. This technique Integrates out numerical difficulties In the lower 
duct region caused by diffusivity gradients and a steep exponential magnetic 
variation. The Integration forces the model predictions to be consistent 
with causality. 

Figure (L-1) shows a representative flux balance at position x in the duct. 

Mass conservation demands that 

REPRODUCffilUTY OP THE 

r>T»T^ r. rn 

2a 

= I in ^ 


X+V2 

r 

x-1/2 


Fy w dx 


Assuming steady state operation over an interval t^ seconds, the amount of 
particulate collected on the lower duct surface is 


/ mass collected \ _ 

' X7F ^ 


x+X/2 


0 x-X/2 


Fy dx dt 


(L-2) 


= w t. 


x+X/2 


-X/2 


r^dx 


With an average axial duct velocity, in the duct, the net mass injected 

avg 


IS 


2a t 


(mass injected) = 


^0 


n^ U w dt dy 
0 


0 0 


(L-3) 


% ^avg ^ ^0 





combining (L-1), (L-3) with (L-2) gives 


^ mass collected ^ 
X/2 


^ mass Injected ^ 
% ^avg 


2a 

0 


(L-4) 


The program predicts n/n^ and since * n u, the mass collected can be 
calculated. The FORTRAN program KDFFl incorporating these changes is listed 
in Appendix M. Simpson's rule is used to perform the transverse flux inte- 
gration. 
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Appendix M - KDFF1 

This FORTRAN program computes densities and half wavelength deposi- 
tion for light particles. Particle deposition Is computed using mass 
conservation arguments; Integration of the entire transverse profile Is 
Involved In predicting deposition. The analysis occurs in Appendix L. 

J CALCULATION OF DUCT CONCENTRATION PROFILE 
J SMALL PARTICLE TURBULENT DIFFUSION CASE 
; BOUNDARY CONDITIONS CONSISTENT WITH CAUSALITY 
I NO BOUNDARY CONDITIONS ON LOWER SURFACE 


dimension b(lOO) i»k(1000) t3n(. 1000) fbO(9) fb2(100) fb2k(100) yy (100) 
dimension ur( 100) i.ariO< 100) 7bb(100i-4) ? a (lOOr 3) y ans ( 2; 100) »ncs( 40) y 
ncf(140) 

double precision ayb 
common ak y bO y m y ammv r nha rm 
wri te(6y5) 

5 Formate THE NUMBER OF ROWS ANIi COLUMNS AF<E") 
read(5yl0)irowf icol 
10 format<2i5) 

write(6yl5)i row y icol 

15 format<2i5y/y ■ INPUT TOyAMMOy WOyXByDLXMS ARE* > 
read ( 5 y 20 ) tsj y ammv y wv y ;cb y dl5<ms 
20 format <fl5.S) 

w r i te ( 6 y 25 ) t d y xb y ammv y w v y dl xms 
25 format(" THE INPUTS TGyXBy AMMVy WVyDLXMS ARE " y/6f 10 . 4 ) 
nhsrm»9.01 
nu=0 

no=ieol/25 
ncol 1=0 

:i. f ( wv ♦ eo ♦ 3 1 . o r . w v ♦ eo ♦ 1 2 . ) do to 30 

ak==15.71 

;<mad~ . 439 

bO(l)™» 080356 

bO < 2 ) =■- . 0027 

b0<3)=!. 00392 

b0(4):=»00151 

b0(5)=»001 

b0(6)-»0000931 

bO.C 7) =-.0007905 

b0<8>™. 0005376 

bO (9 )~. 000352 

do to 40 

30 if (wv.ei3.12. )do to 35 
ak'==9 .975 

>cmad~ .3175 
bOd ) = . 099763 

31 b0(2)-.001S3 
b0(3)«.00457 
b0(4)“. 0005366 


% -Ai 


%• 
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I 


*■ 


b0<S)«»001S 
b0(6)«. 000678 
b0<7)«. 00126 
bO(0)*.OOOS 

bO ( 9 > « . 00095 rEPRODUCBILITY OF THE 

to 40 0RIGINAT> PAGE IS POOR 

35 ak«6*650 

xiri3sJ». 40005 
b0<l)«»08993 

38 b0(2)»*001S6 

b0(3)=»025 
b0(4>«» 000267 
bO(S)«. 00123 
b0<6)=»0000583 
b0(7)»»00364 
b0(8)-. 000267 
b0(9)«. 0001064 
40 dlx a l*/(icol) 

ww»12*7 
wa* 127 

ndl f »wv/2 ♦ 'M ♦ 01/dlK/xb 

kcsaO 

kc«0 

nendaicol-no 
do 43 Jano»nend»no 

42 kcakc+1 

if < <xb>KJ)KdlK“Wv/4H<*01) <kcN<;jb)Kdl;;) > do to 42 

kcsakcs+1 

ncs<kcs)=kc 

ncf(kcs)a kc+ndlf 

43 continue 
kcsendakcs 

write<6i»44)kcserid» (ncs( J) yncf < J) fJ~l okcsend) 

44 formate' kesend k ncs y and nef are ' i3/1.4i5 ) 

kc»l 

dlya ( 1- ♦ 031/ww ) /i row 
do 45 i^iyirow 
y < i ) al ♦ -dlyJl^i 
anOei )=1 f 
45 continue 

do 50 i=lyirow 
call uOam ( y < i ) y am y amk y u ) 
b2 ( i ) aam 
b2k < i ) a amk 
Ur < i ) ai,j/;.'b 
50 continue 

c INPUT FINISHED— SET UP FINITE ELEMENT EQUATIONS 
c KENT LOVES LANA 

do 151 ncolalyicol 
xposaxmad + dl;;:K>!b>Kncol 
dyadly 

do 100 nrowalyirow 

call dt2a < xpos y dt y y < n row ) y dtmax ) 

tidadt/w:K;K2 
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a3"l/u*’<'2/dyn0. 

c=0. 

i f ( y ( n row ) » < 2 1 , » 875 ) c~- “ 1 , 
if <y <nrow) ♦ it* ♦ 125')c=l ♦ 

strdt « c!f.dtiii3.K/w3K!lc2H{t<3Hcl2.7/<1.5875“.031)/2./dw 
u < nrow 1 2 ) =t;3JKt id5ic2/dyH(J((2+ < 1 ♦ f tdH<b2 < n row ) /w ) /dy+ta>Ku r ( nrow ) / 
\c(rirow> 1 dlxftg*b2k(nrow) 

aCnrowr 1 )=-tdJKtid/dy3KNc2-( 1 ♦ +ta>lcb2(rirow)/w)/dy-:^pdt 
a ( n row r 3 ) --t£j-1« t i d/dy5KH<2+£i rdt 
b ( nrow ) r < n row > /dlx^anO < n row ) 

iOO continue 

a ( i row» 1 >-a < i rowy 1 ) fiardt-dtmaxJK ( ♦ 031/ww+dy ) / ♦ 12S./dy>Ktd *d3 
a(irowy2) »3< irowy2)+dtiT(3X'-K< *031/ww+dy) / ♦ 125/dy?i{td *d3 
3(1 y 1 )=0» 
a( irowy3)«0* 

c SINGLE COLUMN OF ELEMENTS NOW READY FOR INVERSION 
139 call invrse(ayby irow> 

do 145 JJ^lyirow 
anO(JJ) ~ b(JJ) 
net) X CO 1/4 

i F ( ncol ♦ ne . nco ♦ and ♦ ncol ♦ ne ♦ 2^nco * and . ncol ♦ no ♦ 3)Knco ♦ and ♦ ncol 
.ne. (4*nco-l)) go to 145 
if (ncoll *ea*ncol )ao to 903 
nu=nu-f 1 

903 ncoll»ncol 

bb ( J J y nu ) =an0 < J J ) 

1<^5 continue 

ridum=l 

if ( ncol ♦ ne » nc3 (kc )♦ and ♦ ncol ♦ ne ♦ ncs ( kef 1 ) ) do to 1S5 
if (ncol ♦de.ncs<kc+l ) ) nduiTi=2 
do 180 i-lyirow 
ans ( nduiTi y i ) -anO ( i ) 

130 continue 

135 X f (nef ( kc ) ♦ eo * ncol )cal 1 f luxxn( ans y anO y x row y ur y dly y dif a ; 

X f ( nef ( kc ) * ne ♦ ncol ) do to 1 5 1 
if ( ncol ♦ It . ncs (kc+1 ) ) do to 170 
do 137 i-lyxrow 
ans ( 1 y X ) - ans ( 2 y ;i. ). 

J.37 continue 


L70 kc^:=kc-fl 

an ncol ) -100 ♦ !{{dif a;«Kb/ . 9/4 * 572 
15 L continue 

c DENSITIES ARE NOW KNOWN - OUTPUT FOLLOWS 
noo-ico.'l/4 

w r X te ( 6 y 90 ) ( i y i -noo y i co 1 y noo ) 

-0 format (" THE DENSITIES ARE AS FOLLOWS 

do 103 i-lyirowy3 
w r i t e ( 6 y 1 07 ) y < i ) y ( bb ( i y J ) y J- 1 y 4 ) 

107 format(f5*3y 4f8*4) 

108 continue 

d o 190 J - 1 y k c s e n d 
X ( X ) -ncs ( i ) ){<dlx!«xb f wv/4 ♦ )K * 01 
190 continue 

wri te ( 6 y 125 ) ( x ( i ) y an ( i ) y i~l y kesend ) 
125 format < 4el5 » 4 ) 

StOF 

end 


y/ylOx 


i 4 y 3 < 6x y i 4 ) ) 
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s».jb routine? uOam ( yu r am f airiK * u ) 
dimansion bOC9) 
common 3ki*b0?bi»3mmv>nh3rm 
y=yu 

if <yu*dt » ♦5>y=l »-y 

U's=4 ♦ 572* ( 2 > *y ) ** ♦ .1. 1 .11 1 1 1 

ainwO ♦ 

do 15 i«l»nharm 
>='=2**i*3k*yu 
if (p»dt *20* )do to 15 
am»am+ammv*bO< i )**2*i*ex»S'<~p) 

3mk»amk+2* *i**2*ammv*b0< i ) **2*3k/w*exjs' <-p > 

15 continue 

if (yu. sit* *0309/12*7) slo to 30 
u == *242**2*yij** 127/ *000015 
30 return 

end 

subroutine dt2a < k a dt » y ? dtma;c ) 
dj{'=*1524 

i f < K * 1 1 ♦ dx ) dt= ♦ 003625 
xel“2 **dx 

if < it ♦ <3e * dx . and * x * It * xel )dt= * 003625 f * 020915* < x~dx ) /dx 
if (x*se* 1*33985) dt** 00178 
if (x*lt*xel*or*x*sie*l*33985) So to 50 

dt* * 054658- * 1 50851 *x+ ♦ 20 1784*x**2- * 084262*x**3- ♦ 07902*x**4 
. l+.028938*x ’ **5+. 065659*x**6- . 033616****7 
SO dtmax dt 

if <y*lt* *125)dt*dt*y/*125 

if (y.le* .031001/12*7) dt * 0. 

if <y *.st* .!:J75)dt*dt*(i *-y)/* 125 

return 

end 

subrou '.r;e invrse < 3 > b f i »'ow ) 

dimension a< 13S?3) fb<13S) yd< 135) »u< 135^2) i»y( 135) 

double precision s?b 

U(l!*l) * 3(1 i»2) 

u(1^2) * 3(1*3) 

do 30 i ” 2 !' i row 

u(i*l) ™ 3(i?2)-3<i-ly3)*a<is.l)/u(i-lyl) 
u ( i y 2 ) * a ( i y 3 ) 
d ( i ) * 3 ( i y 1 ) /u ( i - 1 y 1 ) 

30 continue 

y ( 1 ) * b ( 1 ) 

do 40 i = 2yirow 

y ( i ) * b ( i ) -d < i ) *y ( i - 1 ) 

40 continue 

b ( i row ) * y ( i row ) /u ( i row y 1 ) 
do 50 J - 2yirowyl 
i = irow-fl-J 

b(i) * <y ( i ) ~u ( i y 2 ) *b( i+1 ) )/u( i y 1 ) 

SO continue 

return 
end 
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sub routine f luKin < ans » anO * i row rurrdlyrdifa) 

dimension ans(2i»100) i»anO< 100) Hjr( 100) 

f l«dly/3*5fcans< i row))Kur( i row) 

f2»dly/3«*3n0< i row)!Kur( i row 

i ro®i row-1 

iroe®i row-2 

do 10 i^lfiroi-2 

fl®fl+ur(i ))f<4*^dly/3»’f{3ris<i ) 

f2®f2 f ur(i )/ic4«H{dly/3Jff3n0< i > 

10 continue 

do 20 i~“2yiroe>2 
f l=f 1+ur ( i ) JK2 fifed lu/3 ♦ iKans ( i ) 
f2-f 2+ur ( i ) )K2 fifed ly/3 f ifeanOC i ) 

20 continue 

difa=fl-f2 

return 

end 
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